A Consistent Methodology to Evaluate Temperature and Heat Wave Future Projections for Cities: A Case Study for Lisbon

: Heat waves are large-scale atmospheric phenomena that may cause heat stress in ecosystems and socio-economic activities. In cities, morbidity and mortality may increase during a heat wave, overloading health and emergency services. In the face of climate change and associated warming, cities need to adapt and mitigate the e ﬀ ects of heat waves. This study suggests a new method to evaluate heat waves’ impacts on cities by considering some aspects of heat waves that are not usually considered in other similar studies. The method devises heat wave quantities that are easy to calculate; it is relevant to assessing their impacts and permits the development of adaptation measures. This study applies the suggested method to quantify various aspects of heat waves in Lisbon for future climate projections considering future mid-term (2046–2065) and long-term (2081–2100) climates under the RCP8.5 greenhouse emission scenario. This is achieved through the analysis of various regional climate simulations performed with the WRF model and an ensemble of EURO-CORDEX models. This allows an estimation of uncertainty and conﬁdence of the projections. To evaluate the climate change properties of heat waves, statistics for future climates are compared to those for a reference recent climate. Simulated temperatures are ﬁrst bias corrected to minimize the model systematic errors relative to observations. The temperature for mid and long-term futures is expected to increase relative to the present by 1.6 ◦ C and 3.6 ◦ C, respectively, with late summer months registering the highest increases. The number of heat wave days per year will increase on average from 10, in the present climate, to 38 and 63 in mid and long-term climates, respectively. Heat wave duration, intensity, average maximum temperature, and accumulated temperature during a heat wave will also increase. Heat waves account for an annual average of accumulated temperature of 358 ◦ C · day in the present climate, while in the mid and long-term, future climates account for 1270 ◦ C · day and 2078 ◦ C · day, respectively. The largest increases are expected to occur from July to October. Extreme intensity and long-duration heat waves with an average maximum temperature of more than 40 ◦ C are expected to occur in the future climates. to evaluate HW and their changes in future climates, HIST temperatures are validated against EOBS. This is performed for NC and BC data for WRF and CORDEX models. In addition, the comparison of WRF with the CORDEX ensemble gives an estimate of the uncertainty associated to the representation of present climate temperature and HW statistics and gives confidence to evaluate their future changes.


Introduction
Heat waves (HW) are extreme weather high-temperature events that may pose a risk to ecosystems and socio-economic activities, in general, and particularly to human health [1]. These phenomena are frequently disregarded since they do not cause as dramatic visual impacts as many other extreme

•
Other HW-related variables are considered, namely, the average daily maximum temperature, the absolute daily maximum temperature during a HW, and the sum of the daily maximum temperatures. These properties bring extra information for impact and mitigation studies.

•
Monthly analysis for each extended summer months (May to October) when HWs have the strongest effects. • Analysis by HW type according to its average daily maximum temperature. • Inter-annual variability of HW. • Estimation of return period for some HW properties.
We believe that, although extensive, these set of parameters used to characterize HWs are complementary and useful for the development of adaptation measures. The data used and the applied methods are detailed in Section 2. Section 3 presents the climate statistical properties of HW and their changes for future climate scenarios, and the discussion of the results and the main conclusions are stated in Section 4.

The WRF Model
In this study, the authors use a number of climate simulations performed with the WRF v3.5 (Weather Research and Forecasting) model. WRF is a numerical weather forecast model developed for operational and research purposes fully described by [24].
The model was adapted to perform climate simulations by Marta-Almeida et al. [25], who also extensively evaluated its performance for present climate conditions over the Iberian Peninsula. WRF was implemented with three domains that interact in two-way nesting mode and with horizontal resolutions of 81, 27, and 9 km, respectively ( Figure 1). All model settings are those used by Marta-Almeida et al. [25].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 35 • Other HW-related variables are considered, namely, the average daily maximum temperature, the absolute daily maximum temperature during a HW, and the sum of the daily maximum temperatures. These properties bring extra information for impact and mitigation studies.
• Monthly analysis for each extended summer months (May to October) when HWs have the strongest effects.
• Analysis by HW type according to its average daily maximum temperature.
• Inter-annual variability of HW.
• Estimation of return period for some HW properties.
We believe that, although extensive, these set of parameters used to characterize HWs are complementary and useful for the development of adaptation measures. The data used and the applied methods are detailed in Section 2. Section 3 presents the climate statistical properties of HW and their changes for future climate scenarios, and the discussion of the results and the main conclusions are stated in Section 4.

The WRF Model
In this study, the authors use a number of climate simulations performed with the WRF v3.5 (Weather Research and Forecasting) model. WRF is a numerical weather forecast model developed for operational and research purposes fully described by [24].
The model was adapted to perform climate simulations by Marta-Almeida et al. [25], who also extensively evaluated its performance for present climate conditions over the Iberian Peninsula. WRF was implemented with three domains that interact in two-way nesting mode and with horizontal resolutions of 81, 27, and 9 km, respectively ( Figure 1). All model settings are those used by Marta-Almeida et al. [25].

Simulated Data
Simulations were performed for the recent reference climate (HIST-986-2005) and for two future climate projections namely for the mid-term (MED-2046-2065) and long-term (LONG-2081-2100) future, respectively. These are the same periods as those considered in the 5th Intergovernmental Panel on Climate Change (IPCC) Assessment Report [26]. Future simulations considered the RCP8.5 (Representative Concentration Pathways) greenhouse gas emission scenario [27]. This is one of the most studied scenarios and represents a radiative forcing of 8.5 Wm −2 by 2100 and a continuous increase thereafter [28] corresponding to a global averaged warming of about 4.5 • C in 2100.

Observational Data
Observed data are used in this study with the purpose of validating and correcting simulated data. The E-OBS Version 17.0 dataset (EOBS hereafter) was developed by the European Climate Assessment and Dataset (ECA&D) and consists of a land-only daily gridded dataset for Europe for the 1950-2006 period [36,37]. Data for the 1986-2005 period with a rotated pole grid with a horizontal resolution of 0.22 • were extracted.
All simulated and EOBS data were interpolated to a regular grid with a horizontal resolution of 0.2 • , which allows the comparison between them. Daily maximum (Tmax) and minimum temperatures (Tmin) were extracted for all simulations and for E-OBS, for the grid-point nearest Lisbon.
We decided to use a gridded dataset, such as E-OBS, since these were developed to validate models. E-OBS represents, as model data, the area average whereas station data are local.

Bias Correction
Simulated data may have systematic errors relative to observations resulting from many factors such as poor vertical and horizontal resolution, sub-grid processes that are often semi-empirically parametrized, numerical methods used to solve the equations, and misrepresentation of a lower boundary, amongst others. Therefore, the use of such data in impact, adaptation, and mitigation studies may misrepresent climate change. One way to minimize these biases is to perform bias correction (BC).
Various authors have applied BC in similar studies [38][39][40][41][42][43]. These studies clearly demonstrate that the application of BC methods improve the performance of the models in climate simulation, particularly of extreme values and indices.
In this study, we applied the quantile mapping method of BC [44][45][46][47]. This method corrects the Probability Distribution Function (PDF) of the data for each quantile separately considering the observations (i.e., EOBS) as reference. The procedure used here is described in detail by Amengual et al. [46] and is fully explained by Viceto et al. [23]. This method considers that future PDFs may be different from that of the present climate. BC is applied to all simulated daily data using the three months centered on the month to be corrected. This sub-sampling ensures that some seasonal signature is retained in the process.
Firstly, HIST data are corrected by projecting raw uncorrected data (NC) onto the EOBS PDF. Corrected data for HIST are validated against EOBS data for HIST. Hereafter, NC and BC are used to refer to uncorrected and bias-corrected data, respectively. Secondly, the data for MED and LONG are corrected. The full methodology of the applied bias correction is presented in [23,46].
The EURO-CORDEX platform (http://www.cordex.org/data-access/bias-adjusted-rcm-data/) has BC data for many models, but the BC method is different for each model. This may pose comparison problems, particularly when considering the ensemble of the models. In this study, the BC method quantile mapping is applied to daily minimum and maximum temperature simulated by WRF and to the EURO-CORDEX models (CORDEX hereafter). BC is usually cross-validated by obtaining the transfer functions for a period and validated for a different independent period. Here, a bootstrapping procedure is performed. Ten years of model temperatures are randomly selected from the HIST period. The transfer functions are obtained for this period and applied to the remaining 10 years of the period. The procedure is repeated 50 times. The distributions of the 50 10-year BC temperatures are compared to that for the 20-year HIST period. The objective is to conclude that BC is not dependent on the period selected.

Heat Waves
HW refers to a period of hot weather that may be accompanied by high humidity. They occur due to stationary or slow moving pressure systems and therefore are large-scale phenomena. They have different expressions locally due to other factors such as small-scale surface features (i.e., elevation, surface type, and surface slope orientation). The properties of materials present in cities may also amplify the expression of HW, which is known as the urban heat island effect. The simulated data used in this study is obtained from RCM models, which consider the physical properties of different land-use lower boundary surfaces (i.e., urban, forest, rural, and many others). However, these models do not specifically consider the urban canopy of cities, and therefore, temperature-simulated data may not fully represent the local intensification of HW due to the heat island effect. Nevertheless, since temperatures and HW are compared between climates, this problem is minimized, particularly if the Urban Heat Island (UHI) does not change between climates.
There is no universal accepted definition of HW. HW may be defined using different criteria. Most studies consider that HWs are defined relative to the climate of a particular region, and therefore, the same atmospheric conditions may represent a HW in different locations or not. In a special report of the IPCC [48], a HW was defined as "a period of abnormally hot weather". Most meteorology institutes use the World Meteorological Organization (WMO) definition of HW as: "A marked unusual hot weather (maximum, minimum and daily average temperatures) over a region persisting at least two consecutive days during the hot period of the year based on local climatological conditions, with thermal conditions recorded above given thresholds." [49,50].
Different definitions of HW may be devised depending on the need to focus on a particular sector of application that requires specific thresholds to take actions. Most HW indices consider an absolute or relative temperature threshold value above which the maximum daily temperature must occur for more than a certain number of consecutive days. The Task Team on the Definition of Extreme Weather and Climate Events (TT-DEWCE), established by the WMO Commission for Climatology, presents a detailed discussion and suggests guidelines on the definition and indices of a HW [49].
Most studies consider only Tmax when assessing HWs. Some studies also consider Tmin in the definition of HW [7,51]. This is relevant, since night time temperatures (i.e., Tmin) may represent some thermal relief during a HW [52,53]. Fenner et al. [17] and Zuo et al. [54] present a thoroughly review on the HW definitions and metrics used in HW research.
The choice of the HW definition depends on the objective of the study. Here, we want to supply useful information for society to adapt to climate change in a particular locality. It is not the objective of the present study to evaluate impacts for the thermoregulation of the human body neither to compare HWs changes across distinct locals or regions. Therefore, the main objective of the methodology presented here is to be used locally at other places, without having to change the HW definition. The HW definition is applied equally to the present climate and to future projections considering the present characteristics of HW for the present climate.
Here, we define a HW partially following Russo et al. [4] that is relevant for the objectives of this study. This HW definition was recently applied in [21,23]. HW is defined as a period of at least three consecutive days where the Tmax is above a certain threshold value (Tmax_c) in each day. This climatological threshold is calculated for each day of the year and is equal to the 90th percentile of Tmax of all days on a 31-day window centered on that day. For example, to obtain Tmax_c for each of the 365 days, all 31 days of the 20 years (i.e., 620 days) are considered. As an example, the Tmax_c value for 16 March is calculated considering all days from 1 March until 30 March for the 20 years.
In this study, HW are identified for the period starting on 1 May and ending on 31 October (extended summer), since these are the warmest months. In addition, all HW are considered if they have days within this period even if they begin or end outside it. All HW are identified using Tmax_c calculated for HIST of the respective model, since the objective of this study is to evaluate HW changes for MED and LONG relative to HIST. Only then may this information be used in impact and adaptation studies.
Furthermore, the intensity, duration, and time of the year of HW should also be considered because these features can be relevant to evaluate their impact [7]. Other derived quantities relevant to quantify HW heat stress are calculated. In this study, we consider the following measures related to HW: Clearly, the HW in August causes more stress.

The Urban Heat Island and Heat Waves
HWs in the present study are defined using the maximum daily temperature that occurs during the daytime, typically a few hours after noon. The UHI effect is mainly a night time warming effect. This is generally acknowledged by some studies that also refer to the cooling effects, which are also called an 'inverse' UHI, that the UHI may represent during the day (i.e., UHI < 0 • C) [58][59][60]. This is particularly true for summer, which is the season considered here to study HWs.
The interaction between HWs and the UHI is complex and may depend on the city and the particular HW [61][62][63][64]. To clarify this issue, we performed a set of simulations with the WRF model for the city of Lisbon and for one of the strongest HWs that occurred in the Europe, namely the 2003 HW. The simulations were performed for a very high horizontal resolution of 333 m to more realistically justify the use of an Urban Canopy Model (UCM). The following simulations were performed:
Appl. Sci. 2020, 10, 1149 The UHI intensity is calculated as the averaged urban temperature minus the averaged rural temperature. Figure 2 shows the UHI intensity, averaged for all urban grid points, for this extreme HW for Lisbon for no_UCM, UCM1 and UCM2. One may conclude that: 1.
The UHI in the Lisbon metropolitan area is essentially a nocturnal phenomenon with its intensity reaching between 1.5 and 2.0 • C during the night. However, in certain areas of Lisbon, the UHI intensity at night could reach values greater than 5 • C (not shown).

2.
During the day, urban near-surface urban air temperatures could be up to 0.5 • C lower than its rural surroundings. 3.
The intensity of the UHI is slightly stronger when the UCM is not used. This means that the gross features of the UHI intensity of Lisbon are mainly determined by the three urban surface types already specified in the land-use data of WRF and not by using an UCM coupled to WRF. In fact, the UCM acts to reduce the temperature and accumulated heat when compared to not using the UCM.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 35 2. During the day, urban near-surface urban air temperatures could be up to 0.5 °C lower than its rural surroundings.
3. The intensity of the UHI is slightly stronger when the UCM is not used. This means that the gross features of the UHI intensity of Lisbon are mainly determined by the three urban surface types already specified in the land-use data of WRF and not by using an UCM coupled to WRF. In fact, the UCM acts to reduce the temperature and accumulated heat when compared to not using the UCM. The evaluation of spatial details of the interaction between HWs and the UHI at very high horizontal resolutions (a few hundred meters) in cities may be accomplished using UCMs coupled to regional models and by performing simulations for specific events of only a few days. It is not the purpose of this study to analyze the fine details of the UHI expression during a HW, for example to identify areas within the city where the UHI effect is stronger. Therefore, we do not anticipate significant differences in the temperature and HW-derived metrics mentioned in Section 2.5 between using or not an UCM coupled to the WRF model. We reinforce that this study is concerned with statistical (i.e., climate) changes of HW metrics between climates. In addition, it would be unfeasible, due to computing resource limitations, to perform climate simulations at high enough horizontal resolutions (of the order of a few hundred meters) to simulate all HWs identified in each 20-year climate simulations used in this study with the WRF model coupled to an UCM. In fact, to the authors' knowledge, no one has performed climate simulations with an UCM coupled to a regional model. In fact, many studies have performed WRF-UCM high-resolution simulations for HW events of only a few days.
Therefore, in this study, we do not use an urban canopy parametrization to obtain the fine detail expression of the UHI.

Return Periods
Return periods (RP) are estimated for some of the HW measures. They are computed by adjusting Weibull distribution, whose parameters have been obtained by the moment methods. The The evaluation of spatial details of the interaction between HWs and the UHI at very high horizontal resolutions (a few hundred meters) in cities may be accomplished using UCMs coupled to regional models and by performing simulations for specific events of only a few days. It is not the purpose of this study to analyze the fine details of the UHI expression during a HW, for example to identify areas within the city where the UHI effect is stronger. Therefore, we do not anticipate significant differences in the temperature and HW-derived metrics mentioned in Section 2.5 between using or not an UCM coupled to the WRF model. We reinforce that this study is concerned with statistical (i.e., climate) changes of HW metrics between climates. In addition, it would be unfeasible, due to computing resource limitations, to perform climate simulations at high enough horizontal resolutions (of the order of a few hundred meters) to simulate all HWs identified in each 20-year climate simulations used in this study with the WRF model coupled to an UCM. In fact, to the authors' knowledge, no one has performed climate simulations with an UCM coupled to a regional model. In fact, many studies have performed WRF-UCM high-resolution simulations for HW events of only a few days.
Therefore, in this study, we do not use an urban canopy parametrization to obtain the fine detail expression of the UHI.

Return Periods
Return periods (RP) are estimated for some of the HW measures. They are computed by adjusting Weibull distribution, whose parameters have been obtained by the moment methods. The Weibull distribution was chosen since it gave the best fit to the data when compared to other theoretical distributions such as Gumbel and Gamma distributions. In addition, quantile-quantile (Q-Q) and cumulative distribution functions (CDFs) graphs are shown. The RP estimates performed for the HW measures mentioned above are calculated for each climate, namely HIST, MED, and LONG. For each of these climates, RP do not assume climate change. For example, a RP estimate of 50 years calculated for MED assumes that the MED climate does not change for 50 years in MED. However, the climate change of RP is evaluated by comparing the RP amongst the different periods.

Validation-Recent climate
Before using simulated temperature to evaluate HW and their changes in future climates, HIST temperatures are validated against EOBS. This is performed for NC and BC data for WRF and CORDEX models. In addition, the comparison of WRF with the CORDEX ensemble gives an estimate of the uncertainty associated to the representation of present climate temperature and HW statistics and gives confidence to evaluate their future changes.

Validation-Recent climate
Before using simulated temperature to evaluate HW and their changes in future climates, HIST temperatures are validated against EOBS. This is performed for NC and BC data for WRF and CORDEX models. In addition, the comparison of WRF with the CORDEX ensemble gives an estimate of the uncertainty associated to the representation of present climate temperature and HW statistics and gives confidence to evaluate their future changes.    All models, namely WRF and CORDEX, present a generalized cold bias before bias correction with an all-model average annual bias of −2.5 • C (minimum and maximum biases are −4.2 • C and −1.2 • C). Biases for Tmin are smaller than those for Tmax (not shown). For monthly biases, one model has a bias of −5.9 • C for June/July. These biases are reported in a number of other studies [65] and are, at least, partially related to cold biases present in the forcing global models caused by large-scale dynamics' deficiencies in these models [32,66]  After BC, all models approach EOBS, as expected, since BC projects uncorrected temperature onto the EOBS distribution. WRF presents a small differences relative to EOBS.
The cold bias mentioned before is now more evident with smaller WRF and CORDEX ensemble average biases present for extreme Tmax. Larger inter-model variability is observed for CORDEX for extreme high temperature. Similar results are observed for Tmin (not shown). After correction, by definition, the distributions of model temperature are the same as for the EOBS for HIST and therefore are not shown.
Cross-validation was performed as described at the end of Section 2.4. The 50 10-year bootstrapping runs have an average bias relative to EOBS of −3.3 • C for NC and 0.0 • C for BC. The bias for HIST NC is also −3.3 • C. The 25th and 75th percentiles across the 50 runs biases are −0.18 • C and +0.17 • C, respectively. Therefore, one can conclude that the period chosen to estimate the transfer functions is independent of the calibration period.

Heat Waves
In this section, all HW are identified for HIST, and the above-referred HW measures are calculated for each HW by applying the previously defined HW identification criteria that considers Tmax_c (i.e., climatological percentile 90 of Tmax for HIST). As said before, the HW of all the analysis is performed for the period starting on 1 May and ending on 31 October (extended summer).
The annual cycle of Tmax_c (i.e., threshold) for EOBS, WRF NC, and WRF BC ( Figure 5) presents biases relative to EOBS' similar to Tmax and Tmin. These are very large for WRF NC for spring but, again, after BC differences relative to EOBS are minimized.
the models in the y-axis.
All models, namely WRF and CORDEX, present a generalized cold bias before bias correction with an all-model average annual bias of −2.5 °C (minimum and maximum biases are −4.2 °C and −1.2 °C). Biases for Tmin are smaller than those for Tmax (not shown). For monthly biases, one model has a bias of −5.9 °C for June/July. These biases are reported in a number of other studies [65] and are, at least, partially related to cold biases present in the forcing global models caused by large-scale dynamics' deficiencies in these models [32,66]  After BC, all models approach EOBS, as expected, since BC projects uncorrected temperature onto the EOBS distribution. WRF presents a small differences relative to EOBS.
The cold bias mentioned before is now more evident with smaller WRF and CORDEX ensemble average biases present for extreme Tmax. Larger inter-model variability is observed for CORDEX for extreme high temperature. Similar results are observed for Tmin (not shown). After correction, by definition, the distributions of model temperature are the same as for the EOBS for HIST and therefore are not shown.
Cross-validation was performed as described at the end of Section 2.4. The 50 10-year bootstrapping runs have an average bias relative to EOBS of −3.3 °C for NC and 0.0 °C for BC. The bias for HIST NC is also −3.3 °C. The 25th and 75th percentiles across the 50 runs biases are −0.18 °C and +0.17 °C, respectively. Therefore, one can conclude that the period chosen to estimate the transfer functions is independent of the calibration period.

Heat Waves
In this section, all HW are identified for HIST, and the above-referred HW measures are calculated for each HW by applying the previously defined HW identification criteria that considers Tmax_c (i.e., climatological percentile 90 of Tmax for HIST). As said before, the HW of all the analysis is performed for the period starting on 1 May and ending on 31 October (extended summer).
The annual cycle of Tmax_c (i.e., threshold) for EOBS, WRF NC, and WRF BC ( Figure 5) presents biases relative to EOBS' similar to Tmax and Tmin. These are very large for WRF NC for spring but, again, after BC differences relative to EOBS are minimized.  All HW are identified in each model (i.e., WRF and CORDEX models) by the application of the respective Tmax_c. Figure 6 shows the (a) NDAYS annual average and (b) Tmax_sum HW average for EOBS, WRF CORDEX average, maximum and minimum (vertical bar) for HIST BC. WRF and CORDEX ensemble average present values very similar to EOBS. However, the CORDEX ensemble shows some inter-model variability. These two HW measures are very important to quantify the impacts of HW, since they represent absolute and accumulated temperature relevant to heat stress. NWAVES and DUR may not reflect the same. For example, a small NWAVE number in a given year may be associated to a large HW DUR and/or INT. The same is true for many HW of small DUR each.
All HW are identified in each model (i.e., WRF and CORDEX models) by the application of the respective Tmax_c. Figure 6 shows the (a) NDAYS annual average and (b) Tmax_sum HW average for EOBS, WRF CORDEX average, maximum and minimum (vertical bar) for HIST BC. WRF and CORDEX ensemble average present values very similar to EOBS. However, the CORDEX ensemble shows some inter-model variability. These two HW measures are very important to quantify the impacts of HW, since they represent absolute and accumulated temperature relevant to heat stress. NWAVES and DUR may not reflect the same. For example, a small NWAVE number in a given year may be associated to a large HW DUR and/or INT. The same is true for many HW of small DUR each. INT quantifies the average HW Tmax above Tmax_c and therefore may misrepresent the HW impact in absolute terms. For example, two HW with the same INT may be observed for different Tmax_c. Models simulated an average of approximatelly 10 NDAYS, with each HW contributing with nearly 150 of the accumulated temperature. Figures S1-S3 of the Supplemental material show the same information as Figure 6 but for the remaining HW measures. They show an average of 2.5 HW per year of 4-5 days duration.
Tmax during a HW is on average nearly 2 • C above Tmax_c (INT). The HW daily temperature range (RF) is 14.5 • C on average. The average Tmax_ave of all HW is slight above 33 • C considering all extended summer months (May to October). In each year, HW contribute, on average, with 350 • C·day of accumulated temperature (for example one HW with 10 days and Tmax of 35 • C in each day).

Temperature
Having successfully validated the models relative to EOBS for temperature and HW characteristics, hereafter, only BC data are considered to evaluate climate change.  Having successfully validated the models relative to EOBS for temperature and HW characteristics, hereafter, only BC data are considered to evaluate climate change. Figure 7 shows monthly differences relative to HIST of Tmax for (a) MED and (b) LONG, for WRF and average, maximum and minimum CORDEX models. For MED, all models simulate an annual average warming of 1.9 °C with the greatest Tmax increases from July to October. For LONG, changes relative to HIST are similar but stronger than those of MED with the models simulating an average warming of 3.7 °C with greatest Tmax increases also for the July-October period of 4.4 °C. Q-Q graphs for MED/LONG, relative to HIST, are shown in Figure 8. For MED, all models simulate an annual average warming of 1.9 • C with the greatest Tmax increases from July to October. For LONG, changes relative to HIST are similar but stronger than those of MED with the models simulating an average warming of 3.7 • C with greatest Tmax increases also for the July-October period of 4.4 • C. Q-Q graphs for MED/LONG, relative to HIST, are shown in Figure 8. Warming in MED relative to HIST is greatest between percentiles 50 and 80, and uncertainty (model dispersion) is larger for extreme low and high temperatures. WRF is very similar to the CORDEX ensemble average. Tmin (not shown) increases by 1.6 • C in MED and 3.2 • C in LONG, with the greatest increases observed from July to October. Together, Tmax and Tmin changes represent an increase of the average daily temperature range, relative to HIST, of 0.2 • C for MED and 0.5 • C for LONG.
Warming simulated by individual models is shown in Figure 9 for Tmax. Inter-model variability is about 0.2 • C and 0.5 • C (standard deviation) for annual and monthly conditions, respectively. Warming simulated by individual models is shown in Figure 9 for Tmax. Inter-model variability is about 0.2 °C and 0.5 °C (standard deviation) for annual and monthly conditions, respectively. Cardoso et al. [32] evaluated temperature changes for Portugal and, while not comparable, they indicate similar warming for Portugal by the end of the century.

Heat Waves
HW are identified for HIST, MED, and LONG for WRF and CORDEX models. This HW analysis is performed only for extended summer. Figure 10 presents the annual average of NWAVE (a) and NDAYS (b) for WRF HIST, MED, and LONG; and the CORDEX average, maximum, and minimum for MED and LONG.  Cardoso et al. [32] evaluated temperature changes for Portugal and, while not comparable, they indicate similar warming for Portugal by the end of the century.

Heat Waves
HW are identified for HIST, MED, and LONG for WRF and CORDEX models. This HW analysis is performed only for extended summer. Figure 10 presents the annual average of NWAVE (a) and NDAYS (b) for WRF HIST, MED, and LONG; and the CORDEX average, maximum, and minimum for MED and LONG. NWAVES trends are positive with values increasing from 2 per year in HIST to 7 and 10 per year in MED and LONG, respectively. More relevant is the increase in NDAYS by a factor of 3-4 and more than 6 in MED and LONG, respectively, representing an average of 30-40 and 60 NDAYS in the future climates. In line with these results, Cardoso et al. [32] refer to an increase in the number of extended summer (MJJAS) days with Tmax above the 90th percentile from 15 days in the recent climate to 80 days by the end of the century.
The remaining HW measures are shown in Figure 11, namely, DUR, INT, RF, Tmax_ave, and Tmax_sum (annual average and HW average). All these quantities are projected to increase in the future, more in LONG than in MED, except for the RF (Figure 11c), which reduces by about 0.5 °C. This represents, for the future, reduced relief during the night of HW days when daytime temperatures are expected to be even higher than in HIST (see Tmax_ave in Figure 11d).
Garcia-Herrera et al. [1] studied the impact of extreme temperatures on mortality for Lisbon and identified a critical temperature of 33.4 °C above which there is an increase of 31.3% in mortality for each 1 °C above the threshold. Here, Tmax_ave for the extended summer is 33.7 °C for HIST and 34.1 °C and 34.9 °C for MED and LONG, respectively. The annual average Tmax_sum (Figure 11f), All these quantities are projected to increase in the future, more in LONG than in MED, except for the RF (Figure 11c), which reduces by about 0.5 • C. This represents, for the future, reduced relief during the night of HW days when daytime temperatures are expected to be even higher than in HIST (see Tmax_ave in Figure 11d).
Garcia-Herrera et al. [1] studied the impact of extreme temperatures on mortality for Lisbon and identified a critical temperature of 33.4 • C above which there is an increase of 31.3% in mortality for each 1 • C above the threshold. Here, Tmax_ave for the extended summer is 33.7 • C for HIST and 34.1 • C and 34.9 • C for MED and LONG, respectively. The annual average Tmax_sum (Figure 11f), which combines the average duration and average Tmax_ave, increases from 358 • C·day to more than 1000 • C·day and 2000 • C·day in MED and LONG, respectively. Recent positive trends of HW duration, intensity, and frequency in Extremadura (Spain) reported by Acero et al. [18] are consistent with the results shown in this study. Figure 12 shows annual total NWAVES, NDAYS, and Tmax_sum for (a) HIST, (b) MED, and (c) LONG.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 21 of 35 duration, intensity, and frequency in Extremadura (Spain) reported by Acero et al. [18] are consistent with the results shown in this study. Figure 12 shows For Figure 12, CORDEX models are not shown, since it is not useful and meaningful to show, for each CORDEX model, the respective annual values of NWAVES/NDAYS, as well as their ensemble average. These are yearly average values but, in particular years, values can reach 100 NDAYS in LONG, which represents much of the extended summer under HW conditions. In line with these results, Cardoso et al. [32] refer to an increase in the number of extended summer (MJJAS) days with Tmax above the 90 th percentile from 15 days in the recent climate to 80 days by the end of the century.
Some extreme values of these HW measures values are shown in the supplementary material ( Figure S4). These, together with the inter-annual variability of NWAVES, NDAYS, and Tmax_sum shown in Figure 11, give a glimpse of what may happen in some years/HW in the future.
Another aspect that may be relevant for impact and adaptation studies is the HW distribution in the seasonal cycle (here, only the extended summer is considered). This is shown in Figure 13 for NDAYS (a) and Tmax_sum (b) for HIST, MED, and LONG. For Figure 12, CORDEX models are not shown, since it is not useful and meaningful to show, for each CORDEX model, the respective annual values of NWAVES/NDAYS, as well as their ensemble average. These are yearly average values but, in particular years, values can reach 100 NDAYS in LONG, which represents much of the extended summer under HW conditions. In line with these results, Cardoso et al. [32] refer to an increase in the number of extended summer (MJJAS) days with Tmax above the 90th percentile from 15 days in the recent climate to 80 days by the end of the century.
Some extreme values of these HW measures values are shown in the supplementary material ( Figure S4). These, together with the inter-annual variability of NWAVES, NDAYS, and Tmax_sum shown in Figure 11, give a glimpse of what may happen in some years/HW in the future.
Another aspect that may be relevant for impact and adaptation studies is the HW distribution in the seasonal cycle (here, only the extended summer is considered). This is shown in Figure 13 for NDAYS (a) and Tmax_sum (b) for HIST, MED, and LONG.
Increases are observed in MED and LONG for all months, but the largest increments are observed for the late extended summer (September and October) for WRF and CORDEX ensemble average. However, note that in HIST, this is not observed. The monthly disaggregation of other HW measures is shown in the supplementary material ( Figure S5). These changes may represent harsher thermal conditions for rural neighboring areas since it occurs at the end of the dry season, when vegetation is driest, and may be relevant for water resources, power consumption, and agriculture, but most importantly for rural fires [11]. These projected changes of extreme temperature coexist with an extension of the dry season and less precipitation in summer and autumn [66,67].
Finally, Figure 14 shows the contribution of HW with different Tmax_ave for (a) NDAYS and (b) Tmax_sum. CORDEX models are not shown. For both measures and for LONG, extremely intense HW (i.e., with Tmax_ave equal or greater than 40 • C) contribute with 8% and 10% to the total, respectively. These extreme HW are not observed in HIST and MED. Appl. Sci. 2020, 10, x FOR PEER REVIEW 23 of 35 Increases are observed in MED and LONG for all months, but the largest increments are observed for the late extended summer (September and October) for WRF and CORDEX ensemble average. However, note that in HIST, this is not observed. The monthly disaggregation of other HW measures is shown in the supplementary material ( Figure S5). These changes may represent harsher thermal conditions for rural neighboring areas since it occurs at the end of the dry season, when vegetation is driest, and may be relevant for water resources, power consumption, and agriculture, but most importantly for rural fires [11]. These projected changes of extreme temperature coexist with an extension of the dry season and less precipitation in summer and autumn [66,67].
Finally, Figure 14 shows intense HW (i.e., with Tmax_ave equal or greater than 40 °C) contribute with 8% and 10% to the total, respectively. These extreme HW are not observed in HIST and MED.

Extreme Heat Waves
In the previous section, average HW measures were considered. Here, extremes of HW measures are compared amongst HIST, MED, and LONG. Figure 15

Extreme Heat Waves
In the previous section, average HW measures were considered. Here, extremes of HW measures are compared amongst HIST, MED, and LONG. Figure 15  There is a marked shift of percentiles toward higher values from HIST to MED and LONG. The exception is for RF, where percentiles decrease in the future (MED and LONG), meaning a reduction of average daily temperature range during HW. This reduces the relief that people may experience during the nights of HW days when the Tmax itself is higher in MED and LONG than in HIST. Some extreme HW are identified for the three climates. Although rare by definition, they shed some light on extremely hot conditions that can be experienced, or what may be a mega HW in the future. The criteria used here to identify extreme HW are the longest DUR, highest Tmax ave, and Tmax_max.
In this study, the statistical analysis of HW was only performed for the warmer extended summer, since these may be relevant to heat stress conditions. However, due to warming, it may be that HW occurring in the future, in the colder months (November to April), may also cause heat stress similar to that experienced in the present climate during summer. Since the HW identification performed by the authors considered all months, a sample of extreme HW simulated in the colder months, together with extreme HW for the extended summer was selected. As mentioned before, the models used in this study may underestimate the Urban Heat Island effect in locally amplifying the magnitude of the HW; therefore, temperatures may be higher than shown. Table 2 shows a summary of some HW measures, namely, DUR, Tmax_ave, and Tmax_max, for some of the extremest HWs, in terms of duration/intensity, in each climate, simulated by WRF, for colder months (a) and for the extended summer (b). This list of HW is not exhaustive. These help to anticipate, for example, that HW occurring in the non-extended summer in the future (MED or LONG) may have characteristics similar to those occurring in HIST in the extended summer. There is a marked shift of percentiles toward higher values from HIST to MED and LONG. The exception is for RF, where percentiles decrease in the future (MED and LONG), meaning a reduction of average daily temperature range during HW. This reduces the relief that people may experience during the nights of HW days when the Tmax itself is higher in MED and LONG than in HIST. Some extreme HW are identified for the three climates. Although rare by definition, they shed some light on extremely hot conditions that can be experienced, or what may be a mega HW in the future. The criteria used here to identify extreme HW are the longest DUR, highest Tmax ave, and Tmax_max.
In this study, the statistical analysis of HW was only performed for the warmer extended summer, since these may be relevant to heat stress conditions. However, due to warming, it may be that HW occurring in the future, in the colder months (November to April), may also cause heat stress similar to that experienced in the present climate during summer. Since the HW identification performed by the authors considered all months, a sample of extreme HW simulated in the colder months, together with extreme HW for the extended summer was selected. As mentioned before, the models used in this study may underestimate the Urban Heat Island effect in locally amplifying the magnitude of the HW; therefore, temperatures may be higher than shown. Table 2 shows a summary of some HW measures, namely, DUR, Tmax_ave, and Tmax_max, for some of the extremest HWs, in terms of duration/intensity, in each climate, simulated by WRF, for colder months (a) and for the extended summer (b). This list of HW is not exhaustive. These help to anticipate, for example, that HW occurring in the non-extended summer in the future (MED or LONG) may have characteristics similar to those occurring in HIST in the extended summer.   In October HW with long duration (i.e., 20 days) and Tmax_ave with more than 34 °C.


In Summer, HW with long duration (i.e., 10 days) and Tmax_ave in excess of 42 °C.
In general, HW with high DUR and Tmax_ave are expected to occur in the future during spring and autumn than those that were typically observed during summer. For the future projections, HW of strong intensity and duration, never observed before, are expected to occur.

The 2003 Heat Wave in Lisbon
To have a term of comparison with the extreme HW referred before, the unprecedently strong HW of 2003 is analyzed for Lisbon. This HW is reviewed by Garcia-Herrera et al. [68] and Some considerations are as follows: • In March and November HW with Tmax_ave greater than 29 • C.

•
In October HW with long duration (i.e., 20 days) and Tmax_ave with more than 34 • C.

•
In Summer, HW with long duration (i.e., 10 days) and Tmax_ave in excess of 42 • C.
In general, HW with high DUR and Tmax_ave are expected to occur in the future during spring and autumn than those that were typically observed during summer. For the future projections, HW of strong intensity and duration, never observed before, are expected to occur.

The 2003 Heat Wave in Lisbon
To have a term of comparison with the extreme HW referred before, the unprecedently strong HW of 2003 is analyzed for Lisbon. This HW is reviewed by Garcia-Herrera et al. [68] and contributed to more than 58% excess deaths in Portugal mainland [69] and 40% in the Lisbon district [70].
The following analysis is done using observations of temperature recorded at various meteorological stations within the Lisbon area, which wsa obtained from the Instituto Português do Mar e da Atmosfera (IPMA) (see Table S1 for the names and coordinates of the stations). Tmax_c from EOBS is used to identify this HW. Accordingly, the HW is identified starting on 29 July and ending on 2 August 2003 (i.e., 5-day duration). However, since temperatures were high in neighboring days, the period of 27 July until 14 August (i.e., 19 days) is also considered. The period or similar has also been considered in the above-mentioned studies of the 2003 HW.
The 2003 HW was indeed extreme and compares well with the extreme HW found in simulated data for HIST (see Table 2 and Figure S6 in the supplementary material). This HW recorded a value of Tmax_sum of 200 • C·day, which corresponds to percentiles 88, 73, and 58 of Tmax_sum of HIST, MED, and LONG (see Figure 15e), respectively, meaning that a HW similar to that of 2003 will be the nearly average HW in LONG in terms of Tmax_sum. Cardoso et al. [32] found that the 2003 HW is in the 80th percentile of historical HW amplitudes and refer that, by 2100, more than half of HW will be similar to the 2003 HW.
As said before, the values shown in Figure S6 and Tables S2 and S3 represent local measured Tmax and therefore include the local amplification of the HW by the heat island effect. This may be partially absent from simulated Tmax, since the WRF model configuration does not include an urban canopy parametrization. This amplification effect will be quantified in a subsequent study where the same model configuration will be coupled to an urban canopy scheme.

Return Periods
Here, estimated RP are calculated and shown in Figure 16 for NDAYS/year (a) and Tmax_sum per HW (b), for HIST, MED, and LONG.
To have a term of comparison with the extreme HW referred before, the unprecedently strong HW of 2003 is analyzed for Lisbon. This HW is reviewed by Garcia-Herrera et al. [68] and contributed to more than 58% excess deaths in Portugal mainland [69] and 40% in the Lisbon district [70].
The following analysis is done using observations of temperature recorded at various meteorological stations within the Lisbon area, which wsa obtained from the Instituto Português do Mar e da Atmosfera (IPMA) (see Table S1 for the names and coordinates of the stations). Tmax_c from EOBS is used to identify this HW. Accordingly, the HW is identified starting on 29 July and ending on 2 August 2003 (i.e., 5-day duration). However, since temperatures were high in neighboring days, the period of 27 July until 14 August (i.e., 19 days) is also considered. The period or similar has also been considered in the above-mentioned studies of the 2003 HW.
The 2003 HW was indeed extreme and compares well with the extreme HW found in simulated data for HIST (see Table 2 and Figure S6 in the supplementary material). This HW recorded a value of Tmax_sum of 200 °C·day, which corresponds to percentiles 88, 73, and 58 of Tmax_sum of HIST, MED, and LONG (see Figure 15e), respectively, meaning that a HW similar to that of 2003 will be the nearly average HW in LONG in terms of Tmax_sum. Cardoso et al. [32] found that the 2003 HW is in the 80 th percentile of historical HW amplitudes and refer that, by 2100, more than half of HW will be similar to the 2003 HW.
As said before, the values shown in Figure S6 and Tables S2 and S3 represent local measured Tmax and therefore include the local amplification of the HW by the heat island effect. This may be partially absent from simulated Tmax, since the WRF model configuration does not include an urban canopy parametrization. This amplification effect will be quantified in a subsequent study where the same model configuration will be coupled to an urban canopy scheme.

Return Periods
Here, estimated RP are calculated and shown in Figure 16  As said before, these RP are estimated not considering climate change from the respective climate thereafter. For example, for 40 NDAYS/year, an RP of nearly 200 years is estimated for HIST, not considering climate change. However, considering climate change, for example, from HIST to LONG, they are estimated to be less than 2 years in LONG. Therefore, these estimates should be considered with caution. In addition, they are calculated based on parametric functions fitted to the 20-year data; therefore, RP much longer than 20 years are extrapolated estimates. As expected, the RP of a specific value for the future projections are much reduced relative to HIST.

Conclusions
Cities are densely populated areas where extreme hot conditions may negatively affect vulnerable people and overload health and emergency services. Lisbon is one of the warmest European cities and HW may negatively impact its population. In the face of climate change and associated warming, extreme temperatures, particularly associated to HW, are expected to increase the burden on services and increase morbidity and mortality. Therefore, measures are needed to adapt to the impact of these changes that are occurring at an ever-increasing pace. Most HW studies usually consider the number of HW and HW days, duration, and intensity. The method proposed here devises other HW parameters that are complementary and relevant to evaluate impacts of HWs-namely, the accumulated temperature during an HW and the relief people may experience at night. We believe that, although extensive, the set of parameters used to characterize HWs is most useful for the development of adaptation measures. In addition, the method performs the analysis for each month, evaluates the inter-annual variability of HWs and their properties, and estimates return periods. All simulated data are bias corrected prior to the analysis. This is sometimes ignored in adaptation studies but is highly important since it minimizes model systematic errors.
The method suggested here is applied for Lisbon in terms of what is meant to be relevant for the implementation of practical adaptation and mitigation measures. The simulated temperature data used here are bias corrected prior to the analysis using a quantile-quantile method. This is often not done in most studies where, when applied, some simpler forms of bias correction is applied. This is done for MED and LONG future climate projections under the RCP8.5 scenario relative to a reference present HIST climate that were simulated by the WRF and a 15-member ensemble of EURO-CORDEX models.
All these models presented a marked cold bias relative to EOBS. On average, the bias across all models was −2.5 °C for annual Tmax average, with the largest biases occurring for summer. In order (b) Figure 16. Return period of (a) NDAYS and (b) HW Tmax_sum, for HIST, MED, and LONG.
As said before, these RP are estimated not considering climate change from the respective climate thereafter. For example, for 40 NDAYS/year, an RP of nearly 200 years is estimated for HIST, not considering climate change. However, considering climate change, for example, from HIST to LONG, they are estimated to be less than 2 years in LONG. Therefore, these estimates should be considered with caution. In addition, they are calculated based on parametric functions fitted to the 20-year data; therefore, RP much longer than 20 years are extrapolated estimates. As expected, the RP of a specific value for the future projections are much reduced relative to HIST.

Conclusions
Cities are densely populated areas where extreme hot conditions may negatively affect vulnerable people and overload health and emergency services. Lisbon is one of the warmest European cities and HW may negatively impact its population. In the face of climate change and associated warming, extreme temperatures, particularly associated to HW, are expected to increase the burden on services and increase morbidity and mortality. Therefore, measures are needed to adapt to the impact of these changes that are occurring at an ever-increasing pace. Most HW studies usually consider the number of HW and HW days, duration, and intensity. The method proposed here devises other HW parameters that are complementary and relevant to evaluate impacts of HWs-namely, the accumulated temperature during an HW and the relief people may experience at night. We believe that, although extensive, the set of parameters used to characterize HWs is most useful for the development of adaptation measures. In addition, the method performs the analysis for each month, evaluates the inter-annual variability of HWs and their properties, and estimates return periods. All simulated data are bias corrected prior to the analysis. This is sometimes ignored in adaptation studies but is highly important since it minimizes model systematic errors.
The method suggested here is applied for Lisbon in terms of what is meant to be relevant for the implementation of practical adaptation and mitigation measures. The simulated temperature data used here are bias corrected prior to the analysis using a quantile-quantile method. This is often not done in most studies where, when applied, some simpler forms of bias correction is applied. This is done for MED and LONG future climate projections under the RCP8.5 scenario relative to a reference present HIST climate that were simulated by the WRF and a 15-member ensemble of EURO-CORDEX models.
All these models presented a marked cold bias relative to EOBS. On average, the bias across all models was −2.5 • C for annual Tmax average, with the largest biases occurring for summer. In order to minimize these systematic errors and build confidence in the projections, the same BC method was applied to all simulated daily Tmax and Tmin.
After BC, differences of Tmax and Tmin relative to EOBS practically disappeared, and the WRF agreed well with the CORDEX ensemble. After BC, HW measures were also compared between models and EOBS. Although being temperature-derived quantities that were not directly BC against EOBS, simulated HW measures compare well with those of EOBS, with the WRF practically matching the CORDEX ensemble average.
For HIST, the average NWAVES per year is 2.5, and the average HW DUR is 4-5 days. The average Tmax_ave of all HW is slight above 33 • C considering all extended summer months. This compares with an average Tmax for all extended summer days of 26.2 • C. The Tmin_ave for HW is 19.0 • C and for all extended summer days, it is is 15.5 • C. The HW daily temperature range (i.e., RF) is 14.8 • C on average and 10.7 • C considering all extended summer days. This difference of 4.1 • C represents less relief during the nights of a HW when compared with the typical extended summer day. In each year, HW contribute, on average, with 358 • C·day of accumulated temperature (for example one HW with 10 days and Tmax of 35 • C in each day). Differences of HW measures amongst each month are negligible, except for Tmax_ave with JJA registering highest values (35-36 • C). These are 4-5 • C above the critical temperature of 33.4 • C for Lisbon, above which there is an increase of 31.3% in mortality for each 1 • C above the threshold [1].
Future projections show an increase, relative to HIST, of annual average Tmax of 1.9 • C and 3.7 • C for MED and LONG, respectively. For Tmin, the values are 1.6 • C and 3.2 • C, respectively; these are slightly below those for Tmax, which translates into a slight increase of average daily temperature range. Both temperatures show larger increases from July to October (2.6 • C and 5.3 • C for Tmax in July, in MED and LONG, respectively).
For HW, the annual average of NDAYS is expected to increase more than 3-fold (38 days) and 6-fold (63 days) in MED and LONG, relative to HIST (10 days), respectively. In individual years, NDAYS in excess of 100 are expected to occur, which represents much of the extended summer (180 days) under HW conditions.
In addition, Tmax_ave is expected to increase, on average, from 33.7 • C in HIST to 34.1 • C in MED and 34.9 • C in LONG but with individual HW reaching Tmax_ave values greater than 40 • C in the future, particularly in LONG. For example, the percentage of HW days with Tmax equal or greater than 37.5 • C, relative to the total HW days is 8%, 15%, and 23% for HIST, MED, and LONG, respectively. Extremely intense HW with Tmax_ave above 40 • C occur only in LONG.
The annual average of Tmax_sum increases from 358 • C·day in HIST to of 1270 • C·day in MED and 2078 • C·day in LONG with some years registering more than 3000 • C·day (i.e., equivalent to 85 NDAYS with Tmax of 35 • C in each HW day). It has been shown in Section 4 that HW similar to the one of 2003, which was considered one the most intense in Lisbon, are expected to be an average HW in LONG climate.
Contrary to HIST, the months from July to October have the highest values of NDAYS and Tmax_sum in MED and LONG, with September and October being the highest of all. This change of extreme heat toward the latter months of the extended summer may be critical for water resources, agriculture, and forest fires, since it coincides with the end of the dry season. As mentioned before, other studies have shown that precipitation is also expected to decrease in summer and autumn, in the future, relative to present climate conditions. All these changes must be considered to implement adaptation and mitigation measures to minimize the effects of heat stress in Lisbon.
RP estimates of HW measures were calculated for each climate. As expected, the RP of a specific value for the future projections are much reduced relative to HIST. These may be useful for various socio-economic sectors and may be applied either considering climate change (i.e., from HIST to MED and LONG) or for minimized warming if other drastic greenhouse emission scenarios are implemented.
As mentioned before, the models used here may not fully account for the amplification of the HW effects by the Urban Heat Island effect. These models have detailed high-resolution land-use specifications that consider, amongst others, sub-classes or urban surface. However, to the authors' knowledge, none of these models used an urban canopy scheme. A series of simulations with WRF coupled to an urban canopy parametrization constitute a subsequent step to be performed by the authors to determine, more realistically, the local amplification of HW in Lisbon, particularly during the nights of HW days.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/10/3/1149/s1, Figure S1-Annual average of (a) NWAVES and (b) DUR for EOBS, WRF CORDEX average, maximum and minimum (vertical bar) for HIST BC, Figure S2- Figure S6-Tmax during the 2003 extended HW (i.e., 19 days) in Lisbon. The 'envelope' shows maximum and minimum Tmax recorded by the nine meteorological stations. The grey period represents the 5-day HW identified using the criteria used in this study. The black line is Tmax_c from EOBS, Table S1-Meteorological stations operated by IPMA and respective ID and coordinates for the Lisbon area,