Spatio-temporal Changes in the Heatwaves and Coldwaves in Spain (1950-2018): Influence of the East Atlantic Pattern

The Iberian Peninsula has a complex orography, which determines an important altitudinal gradient and alternation of valleys and mountains, and periodic cold/warm advections air. In the present investigation the evolution of the characteristics of heatwaves (HWs) and coldwaves (CWs) (number of events, frequency, duration, magnitude, and amplitude) was analyzed. A total of 28 homogeneous-period weather stations (1950-2018), grouped into six regions (cluster). After submitting the meteorological series to a process of homogenization and data quality control, various ET-SCI indices were estimated in order to obtain evolution trends in each climatic region. In all cases, there was an increase, often significant, in the recurrence of HW events (0.3 / 10 yrs) as well as a decrease in CW events (-0.2 / 10 yrs). In addition, the evolution of the above indices and anomalies was correlated with the evolution of the global index of the East Atlantic (EAi).


Introduction
Heatwave (HW) and coldwave (CW) events are some of the most impacting situations of the global climate. These are extensive, and sometimes very intense, advections of hot/cold air that reach extreme temperatures over a long period of time, after exceeding predetermined thresholds for each territory. This type of situation has seemed to evolve periodically during recent decades, linked to the consequences of climate change (Cubasch et al., 2013;Liss et al., 2017;Mora et al., 2017).
In Spain, according to Meteorology Statal Agency (AEMET) studies, there are 59 cold waves and 57 heat waves in the period 1975-2019 (Chazarra et al., 2020). The attributable risk (%) associated with mortality due to cold waves in Spain (2000Spain ( -2009 reaches 26% in regions of the center-south of the peninsula, and between 19 and 21% for heat waves in western provinces (Linares et al., 2017). Extreme temperature events have been observed on a recurring basis in recent years, especially HWs. Their analysis has been given a major boost at the global and regional scales, due to the growing concern about global warming (Huth et al., 2000, Beniston et al., 2007, Chauvin & Denvil, 2007Vautard et al., 2007). There is an extensive literature on the health impact and additional effects of such events (D'Ippoliti et al., 2010;Anderson & Bell, 2011;Gasparrini & Armstrong, 2011;Barnett et al., 2012). At the local scale many cities have also had to support HWs, which are often intensified due to the existence of urban heat islands (Basara et al., 2010, Gabriel & Endlicher, 2011, Li & Bou-Zeid, 2013. In some major European cities, the population, not accustomed to extremely high temperatures, is especially vulnerable to their effects -such as in Par-is in August 2003 (Lemonsu et al., 2015). Until now, numerous indices based on the apparent temperature (Kalkstein & Valimont, 1986;D'Ippoliti, 2010) and percentile-based indices (Meehl & Tebaldi, 2004;Fischer & Schär, 2010) have been used to evaluate this type of episode.
However, the cold events have been less studied in recent years, and this work has mainly focused on the proposal and application of cold indicators of cold and heat waves. Similar indices were suggested by Frich et al. (2002) and Keevallik and Vint (2015) to detect climate change. These works generally show tendencies towards a decrease in the occurrence of cold indicators, although there is no common position in the scientific community.
Furthermore, extreme events (HWs or CWs) do not experience the same space-time evolution throughout the planet, as has been addressed in recent decades (Karl et al. 1996;Karl & Knight, 1997;Karl & Easterling, 1999;Meehl & Tebaldi, 2004). Liu et al. (2018) showed that the annual number of frost days did not follow a sufficiently clear pattern, globally, during the period 1982-2012, with decreases in northern Siberia, the Tibet Plateau, and northwestern North America, while a certain increase was observed in some European regions. In fact, changes in atmospheric circulation patterns are altering in an uneven way around the planet the characteristics of extreme thermal events, concerning their frequency, magnitude, and duration (Barriopedro et al., 2006;Dong et al., 2008;Garcia-Herrera et al., 2010). The connections between extreme temperature indices and large-scale circulation patterns, such as the North Atlantic Oscillation (NAO) and the East Atlantic (EA) and East Atlantic/ West Russia (EA-WR) patterns, have been analyzed by various authors, with different results according to the areas studied in Europe (El Kenawy et al., 2011;Burić et al., 2014;Doderović & Burić, 2015;Arsenovic et al., 2015;Ciarlo & Aquilina, 2016;Trbić et al., 2017;Milosevic et al., 2017;Burić et al., 2019).
The objective of the present research is to determine the heat waves (HWs) and cold waves (CWs) trends in the southwest of the European continent (Iberian Peninsula) from the analysis of extensive and reliable meteorological series. For this purpose, indices belonging to a set of extreme climatic indicators proposed by Nairn and Fawcett (2013) and Perkins and Alexander (2013).
In addition, another objective is to analyze the connection between the East Atlantic index (EAi), as an influential large-scale circulation pattern, and the ET-SCI (Expert Team on Sector-Specific Climate Indices) regarding HWs and CWs. The EAi is an index that marks the latitudinal mobility of the subtropical high pressure belt in relation to Western Europe. This is an index similar to the NAO index, although its center of action is more the southeast of the Atlantic Ocean and its dipole in the pressure field is between Iceland and the United Kingdom (Barnston and Livezey, 1987;Murphy and Washington, 2001). The EAi, in its positive phase, reflects positive anomalies in the subtropical North Atlantic during winter, in surface atmospheric pressure (SLP) and in the Geopotential 500 hPa (Z500). Hence, the influence of this index on temperatures in Western Europe has been corroborated (Moore et al., 2012). As indicated by other authors, the EAi is more reliable than the NAO index in the correlation with peninsular temperatures (Lorenzo et al., 2008).

Data and methods
The study area refers to the Iberian Peninsula ( Figure  1), a transition zone between both extratropical and subtropical influences and Atlantic and Mediterranean climates. The climatic variability ranges from temperate climates with regular Atlantic rainfall (2000 mm/year) to the southeast semi-arid climates with less than 200 mm rainfall, with a significant thermal gradient of up to 20ºC in the annual average temperatures (Gómez-Zotano et al., 2015).

Observation data
In the thermal analysis, the daily maximum and minimum temperatures  of 28 meteorological stations belonging to the State Meteorological Agency of Spain (AEMET) and the National Meteorological Service of France (Météo-France) were used (Table 1). In addition, the pattern teleconnection (PT) of the EAi was estimated, using monthly data from 1950 to 2018 (https://www.cpc.ncep.noaa.gov/data/teledoc/ ea.shtml).

Data quality control (DQC)
As a previous step of information quality control a process of homogenization and subsequent reconstruction of the meteorological series, including data gaps, was performed using the homogen function (Eq. 1) of the Climatol 3.1.1 package (https://CRAN.R-project.org/package=climatol) (Guijarro, 2011;Guijarro 2018). This function includes missing data filling and detection and correction of outliers and shifts in the mean of the series. Input data were provided in two text files, one with the daily temperature values and another with the station coordinates. The base data fields were the acronym of the analyzed variable, the initial and final years of the data series as set in the parameters of the following call: Homogen ('variable acronym', initial year, final year, snht, snht2) where snht is a inhomogeneities by moving windows and snht2 the general inhomogeneity of the series Climatol applies a cluster analysis in its initial checks of the data, but the number of clusters is automatic, according to the final Root Mean Squared Error (RMSE) and the standard normal homogeneity test (SNHT). The RMSEs are calculated by comparing the estimated and the observed data in each series and the dendrogram of clusters thus obtained is based on their correlation coefficients. Six groups of stations have been produced at a dissimilarity level of  1.20 ( Table 2). The homogeneous series from stations sharing a common climate (temporal evolution) were clustered, resulting in six different regions.

Trend methods
In order to determine trend patterns in Spain, standardized indices belonging to the Expert Team on Sector-specific Climate Indices (ET-SCI), and subject to an annual analysis, that refer to the magnitude, amplitude, frequency, and duration of the heat waves (HWs) and cold waves (CWs) were adopted here ( Table 3). For their calculation, the ClimPACT2.0 software program (GUI R Studio software) of the Pacific Climate Impact Consortium (PCIC) of the University of Victoria (Australia) (Alexander & Herold, 2016) was applied to the datasets of daily temperatures https://www.wmo.int/pages/prog/wcp/ccl/opace/opace4/ ET-SCI-4-1.php . This is an R software package that calculates the ET-SCI HW and CW indices, as well as additional climate extreme indices, from data stored in text or netCDF files. The algorithms included in ClimPACT2 to estimate the HW characteristics are based on Perkins and Al-exander (2013), hereafter PA13, with some slight modifications to the EHF (Perkins, 2015). The Excess Heat Factor (EHF) is an intensity measure that categorises heatwaves by their severity. The calculation of the EHF (ºC2) is based on a three-day averaged daily Tmean, in relation to the 95th percentile of long-term average temperatures, and the recent (prior 30-day) temperatures, for a particular location (Nairn and Fawcett, 2014). The EHF (Eq. 4) consists of a combination (Eq. 2 and Eq. 3) of two excess heat indices (EHI): where TMi represents the average daily temperature for day i and TM95 is the 95th percentile of TM, which is also calculated within a user-specified base period, over the calendar year and using a 15-day running window. TM = (TX + TN)/2, where TX is maximum daily temperature and TN is minimum daily temperature, in the meteorological summer (JJA). For each index, a map was obtained with the temporal trend by seasons , the temporal evolution by cluster, and an evolution of the standardized anomalies.
The calculation of CWs was performed through the ECF (Excess Cold Factor), developed by Nairn and Fawcett (2013) (Eq. 7). The calculation of the ECF EHF (ºC2) is based on a three-day averaged daily Tmean, in relation to the 5th percentile of long-term average temperatures, and the recent (prior 30-day) temperatures, for a particular location The ECF is a combination (Eq. 5 and Eq. 6) of two cold indices (ECIaccl.) and (ECIsig.), which represent the degree of acclimatization and the climatological significance of cold, respectively. The ET-SCI CW indices were derived from studies on human responses to extremely cold climates. According to these, three consecutive very cold days are needed to significantly increase the mortality rate above its previous value (Wang et al., 2016, Piticar et al., 2018. ECI (sig.)= [((Tmi +(Tmi-1) +(Tmi-2)))⁄3]-Tm05i (6) where Tmi is the average daily temperature of day i and Tm05i is the 5th percentile of Tm estimated for each day (i) of the season considered (from November to March). The ECF (Eq. 5) is defined as: The ECIsig. index measures the degree of excess cold, while the ECIaccl. measures cold stress. Negative ECF values indicate CW conditions, and a period of at least three consecutive days with negative ECF values defines a CW event.
Finally, for the analysis of the EA index, and its correlation with the temporal evolution of the indices, the Kendall concordance coefficient (Kendall, 1938) was adopted, using a significance level (p-value) of 0.99. The monthly data of the EA index were obtained from https://www.cpc.ncep.noaa.gov/.

Results and discussion
Results from the homogenization of temperature data The poor quality of the meteorological series -with erroneous data, missing data, changes in the location of the weather checkpoints, and changes in meteorological sensors-generates meteorological series with a large number of breaks (Acquaotta & Fratianni, 2014). The stations chosen here comply with the consensus of the scientific community (Klein Tank et al., 2002), with maximum and minimum temperature values for at least 80% of the analyzed daily data ( Table 4).
The exploratory analysis of the meteorological series with Climatol showed relatively high SNHT values in some cases. This forced the choice of more aggressive thresholds, of snht1 (140) and snht2 (800) for maximum temperatures, and snht1 (150) and snht2 (1100) for minimum temperatures in the "homogen" function (see Material and methods). In the analysis of the maximum temperatures, some of the series analyzed, such as those of Tarifa (22) and Zamora (26), have had up to nine ruptures since 1950, which led to an intense homogenization process (Figure 2). However, the majority of the stations hardly registered ruptures, with inhomogeneity values from 0 to 1 in 25 of the 28 study stations -as exemplified by Zaragoza (27) and Tarbes (21), whose series hardly contained aberrant anomalies (Figure 3).
After the homogenization process, a drastic reduction in the SNHT values was observed, with an average percentage of original data (POD) of 87.5% (considered quite acceptable). For the vast majority of stations there was hardly any change, resulting in POD values between 95 and 100%. However, some had very low percentages of POD, among which Navacerrada (14), Tarifa (22), and Zamora (26) (20-62%) stand out (Table 4).
The meteorological series of minimum daily temperatures showed fewer ruptures or inhomogeneities: only the stations of Tarifa (22), which had four ruptures, and Ciudad Real (5), with two, exhibited more than one inhomogeneity. The number of inhomogeneities was 0 or 1 for the majority of observatories (Table 4). Therefore, the homogenization process (Climatol) was less complex for the minimum temperatures, since the software had to deal with twice as many total breaks for the maximum temperatures (25) compared to the minimum temperatures (12).
The number of breaks in the minimum temperatures was lower than for the maximum values. This could be due to the greater influence of the changes in the meteorological booths on the maximum values, since their conditions are especially sensitive to the presence of solar radiation. For the maximum temperatures the breaks were distributed homogeneously during the analysis period , although they were drastically reduced after 2005 onwards (Figure 4).

Spatio-temporal trends of the ET-SCI HW and CW indices
Spatio-temporal heatwaves (HWs) variations and trends are described based on the indices number (EHF_HWN), frequency (EHF_HWF), duration (EHF_HWD), magnitude (EHF_HWM), and amplitude (EHF_HWA). The number of heat waves events (EHF_HWN) significantly increased, with 0.3 events/10 yrs in numerous areas, especially in the northern mediterranean coast (C1) and western atlantic area (C2) ( Table 5). Despite the rise in HW events, the average temperature (EFH_HWM) does not seem to have changed excessively during recent decades (0.3ºC2 / 10 yrs). In the northern mediterranean coast (C1) and western atlantic area (C2) west there was a somewhat more significant increase (0.4ºC2/10 yrs).
The number of days of heatwaves per year (EHF_ HWF) has significantly increased since 1950, at an overall rate of 1.9 days / 10 yrs. This has been especially relevant in the observatories C1 (northern Mediterranean coast) and C2 (western Atlantic area), with a rise of up 2.8 to 3.2 days / 10 yrs (Table 5).
The EHF_HWD tended to increase, becoming longer (0.7 days / 10 yrs), although without a statistically significant change. The increase was greatest in the northern Mediterranean (C1) and southern mediterranean area (C4) regions, again on the coast, and also for some observatories of the peninsular center (Madrid), with values of 1.0 to 1.1 days / 10 yrs (Figure 5).
Finally, the amplitude, or maximum peak of absolute temperature in the hottest HW (EHF_HWA), has increased during the last few decades, a statistically significant way (1.5ºC2 / 10 yrs). The increase was significant on the northern Atlantic area (C3) and southern mediterranean area (C4), and the observatories of Santiago de Compostela and Jérez de la Frontera, with an increase of 2.9ºC2 / 10 yrs.
The parameters that refer to the characteristics of coldwaves (CWs) have undergone significant changes during the last eight decades. The number of events (ECF_HWN), since 1950, has decreased significantly throughout the Iberian Peninsula (-0.2 events/10 yrs), especially in the C1 cluster (northern Mediterranean coast) ( Figure 6).
However, the magnitude of the CWs (ECF_HWM) did not undergo many changes, the average temperature of these events showing only a slight rise (0.2ºC2 / 10 yrs) in recent decades. In addition, there were territorial differences, between more marked ascents (C3 in the northern atlantic area, with 0.5ºC2 / 10 yrs) and descents (-0.1ºC2 / 10 yrs, in the C2 western Atlantic  area), the latter being recorded especially in the observatories of the peninsular southwest (Tarifa and Jérez) (Figure 6). The frequency, or number of days of CWs per year (ECF_HWF), showed a statistically significant decrease of -1.6 days / 10 yrs, being generally uniform throughout the study area, although with more marked drops in C1 and C4 (northern Mediterranean coast and southern Mediterranean coast).
The duration of the CW events (ECF_HWD) has also altered during recent decades, with the same evolution pattern as the ECF_HWF. They underwent a general shortening (an average of -0.6 days / 10 yrs), which was more important in the southern mediterranean coast (C4) and the western-central interior region (C5). Specifically, the greatest rates of shortening (ECF_HWD < -0.8) were in Ciudad Real, Barcelona-Fabra, Madrid-Retiro, Burgos and Murcia-San Javier. In addition, the amplitude of the CWs (ECF_HWA) has shown an important increase during the last few decades (2.4ºC2 / 10 yrs), especially in some observatories located in C3 (Bilbao, San Sebastian, and Tarbes), C5 (Burgos), and C6 (Zaragoza) (Figure 6).
In summary, the HW indices show a greater increase in the frequency, duration, and intensity of HW events in the northern Mediterranean (C1) and western Atlantic area (C2), specifically at observatories where there is a maritime influence.
For his part, the ET-SCI CW indices showed a homogeneous behavior, with a trend towards less frequent and less intense CWs; although, it was in the northern Atlantic region (C3), where the greatest increases in amplitude and magnitude occurred, while the western Atlantic region (C2) was the one that experienced the smallest changes (Table 5).

Relationships between the EAi and the ET-SCI HW and CW indices
The results of the present study demonstrate that the EAi has a significant relationship with the temperature records in this territory, which is consistent with the correlations found in other areas of southern Europe (Arsenovic et al., 2015;Buric et al., 2018). The EAi, from 1950 to 2018, showed a temporary upward trend during the winter months (0.1/10 yrs), 1977 being quite significant as a year of rupture, according to the SNHT homogenization test. In addition, during the summer months, the trend was similar (0.1/10 yrs) and also had a statistical consistency at a level of significance of 0.05. In recent years, there has been a winter migration of the jet stream from the North Atlantic to more northern regions (200-300 hPa) (Comas-Bru et al., 2016), which is reflected by above-average surface temperatures in southwest Europe (Wallace and Gutzler, 1981).
Since 1950 The analysis of the relationships between EAi and HWs indicator, by observatory, revealed a higher correlation for those located on the Mediterranean and atlantic coast (Reus, Barcelona-Fabra, Tortosa, Valencia and San Javier) or on the southwestern Atlantic coast (Tarifa and Jerez), and a very low correlation for the northern observatories ( Figure 8). With regard to the CWs, the spatial distribution pattern of Kendall's τ coefficients in relation to the EAi is quite similar, although with slight nuances. The observatories with the highest correlation coefficients contin-     ued to be located exclusively on the Mediterranean coast (from Barcelona, Alicante and Valencia), including pre-coastal areas in the peninsular Southeast (e.g. Murcia-Alcantarilla). By contrast, the worst correlations were obtained in the northern sub-plateau (Castilla and León) and the Ebro basin (Huesca, Zaragoza, and Logroño). Generally, in all the regions derived from the cluster analysis there was a significant correlation (with a confidence level of 0.05) between the evolution of the EAi and that of any ET-SCI index ( Table 6).
The highest average correlations were found for the number of events, frequency (number of days), and duration, with values up to -0.86 in C4 (ECF_HWF) ( Table 6). Considering the relationships between the EAi and the ET-SCI indices by region, it is worth highlighting the highest correlation (-0.82) for ECF_ HWN in C2 and C4 and for EHF_HWF in C4.

Discussion and Conclusion
The homogenization, quality control, and completion of the daily temperature series for the Iberian Peninsula during the period 1950-2018 improved substantially the quality of these datasets and thus our ability to detect the characteristics and trends of the HWs and CWs that have occurred recently in this territory.
The spatio-temporal analysis of the ET-SCI HW and CW indices led to the identification of different patterns of evolution according to the six regions previously defined by clustering. The HWs, although not showing important variations in their average temperature (which, in fact, decreased slightly in observatories in the north of the peninsular), experienced an increase in the peaks of their absolute temperature values, especially in western Atlantic observatories (Santiago de Compostela and Tarifa) and in the peninsular center (Madrid-Retiro). Previously, a similar trend of the HW ratio (0.56 per decade) was shown by Labajo et al. (2014) in the central area of the Iberian Peninsula during the period 1961-2010, and, on a larger scale, by Lhotka and Kyselý (2015) in Central Europe since 1990. Nevertheless, during this century the HWs have tended to be more frequent, longer, and more intense in a large area of the Iberian Peninsula, as has occurred in the rest of southern Europe (Schär et al. 2004;Clark et al. 2006;Alexander, 2010;García Herrera et al., 2010).
Compared with studies conducted in countries with a European continental climate, we have observed here a similar pattern for the HWs (more frequent, longer, and more intense). The seasonal HW trends obtained in the west-central area of the Iberian Peninsula are comparable with those found by Spinoni et al. (2015) in the Carpathian Region, and more particularly in the Hungarian Plain and Southern Romania, where the most severe HWs occur in summer. In addition, most of the HW trends defined for the Iberian Peninsula in this study are quite consistent with the patterns found by Kuglitsch et al. (2010) in the eastern Mediterranean region (e.g. Balkans, Greece, Turkey). They are mainly in accordance with statistically significant increases since the 1960s in the HW number (HWN95), HW length (HWL95), and HW intensity (HWI95). In contrast, the average temperature estimated for this type of event does not seem to have changed as much as in the eastern Mediterranean.
During the last few decades the number of CWs has decreased significantly in the Mediterranean region (C1, C4), including the Ebro basin (C6), and on the western Atlantic coast (C2). Regarding the magnitude, two trends were observed: 1) a little variation in the average temperature of EHF_HWM, and 2) an increase in the absolute minimum temperatures recorded (ECF_HWA) in the northern peninsula (C3) and in observatories of the interior of the Ebro basin (C6) and Barcelona (C1).
Furthermore, a significant relationship between the EAi and the ET-SCI HW and CW indices and their standardized anomalies was corroborated in relation to the number, duration, and frequency of each type of event. The closest relationships with the evolution of the EA pattern, based on correlation coefficients from -0.86 to -0.82, were obtained in the Southwestern mediterranean areas and coastal western atlantic area and for the parameters ECF_HWF and ECF_HWN. Consequently, from all of the above, it is inferred that the temporal evolution of the EAi could be a good indicator for the analysis of future trends in the number of HW/ CW events and their frequency (number of days), associated with climate changes, especially on the Mediterranean coast and western atlantic area of the Iberian Peninsula.
Coldwaves have become less frequent but exhibited variable intensity patterns, according to the Peninsula region, with a significant increase in the average CW temperature in the northern zone and a slight decrease in the southeastern zone under the Atlantic influence. It is striking that in the period 1972-1982 the variations in the HW indices were the slightest of the last 70 years. The same finding was made by Lhotka and Kyselý (2015) when they analyzed the daily temperature series throughout this period in numerous observatories of Central Europe. This behavior was also evidenced in various countries of Central Europe, such as Poland (Tomczyk 2015), the Czech Republic (Kyselý 2010), and Germany (Tomczyk & Sulikowska, 2018).
Many of them correspond to cold periods identified by Labajo et al. (2014) for different meteorological observatories in the center of the Iberian Peninsula during the period 1961-2010. In general, the progressive decline in CWs throughout this period coincides with such a decline in eastern Europe -in both the north, such as Poland (Wibig, 2018;Tomczyk and Bednorz, 2014), and the south, such as areas of the Carpathian region (Spinoni et al., 2015). In winter 2009/2010, a series of extreme cold spells hit northern and western Europe (Cattiaux et al., 2010).
The influence of the EAi on the variability of winter weather in Europe has been widely demonstrated in recent years, linked to the geographical positions of the NAO dipole. When the EA is positive, the action centers migrate to the northeast, while, during the opposite sign phases, they move to the southwest. The EA is, therefore, an excellent indicator of the latitudinal position of the jet stream. With a positive EAi, there is an intense flow of zonal winds from the west over the Iberian Peninsula, with an increase in thermal values in southwest Europe (Mikhailova & Yurovsky, 2016). This corroborates the results obtained here, showing an increase in the frequency and intensity of HWs in summer, and a decrease in CWs in winter in the study area.