Attributing trends in extremely hot days to changes in atmospheric dynamics

. This paper presents a method for attributing regional trends in the frequency of extremely hot days (EHDs) to changes in the frequency of the atmospheric patterns that characterize such extraordinary events. The study is applied to mainland Spain and the Balearic Islands for the extended summers of the period 1958–2008, where signiﬁcant and positive trends in maximum temperature ( T x ) have been reported during the second half of the past century. First, the study area was split into eight regions attending to their different temporal variability of the daily T x series obtained from the Spain02 gridded data set using a clustering procedure. Second, the large-scale atmospheric situations causing EHDs are deﬁned by circulation types (CTs). The obtainment of the CTs differs from the majority of CT classiﬁcations proposed in the literature. It is based on regional series and on a previous characterization of the main atmospheric situations obtained using only some days classiﬁed as extremes in the different regions. Three different atmospheric ﬁelds (SLP, T850, and Z500) from ECMWF


Introduction
The climate suffers changes at different time scales driven by several external and internal factors. Human-induced changes in greenhouse gases, land use, etc., have been especially prominent in the last centuries, modifying the energy balance and therefore inducing climate changes (Stocker et al., 2013). The attribution of recent climate change to each factor is an attempt to ascertain the causes for recent changes observed in the Earth's climate and quantify their relative relevance. However, although the main factors perturbing the climate at a global scale have been extensively characterized (Huybers and Curry, 2006;Swingedouw et al., 2011), fewer exercises focusing on regional scales are available (Stott, 2003).
Beyond the average state, the footprint of climate change is manifested through shifts in extreme weather. In recent years there has been an increasing interest in quantifying the role of human and other external influences on climate in specific weather events (Stott, 2003;Zwiers et al., 2011). Therefore, trying to attribute extreme events to climate change at regional scales presents particular challenges for science.
In the last decades, Europe has experienced a prominent increase in the occurrence of extremely warm episodes, especially during summertime (Frich et al., 2002;Klein Tank and Können, 2003;Alexander et al., 2006). The impacts of such events on human health are important, observing high rates of mortality when such extremes occur. Some examples are the summers of 2003 (Trigo et al., 2005) and 2010 (Dole et al., 2011), when persistent episodes of high maximum temperatures in western and eastern Europe took place. Several studies show that days exceeding the 95th percentile of the daily maximum temperature (T x ) series, usually called extremely hot days (EHDs; García-Herrera et al., 2001), cause also an increase in mortality, especially among the elderly and people with cardiovascular diseases (Díaz-Jiménez et al., 2005). Many works have tried to investigate the causes of these trends, concluding that these are largely influenced by anthropogenic activity (Stott, 2003;Zwiers et al., 2011), with summers like that of 2003 expected to become more frequent under several climate change scenarios (Beniston, 2004). However, extreme events are also driven by unpredictable internal variability. Indeed, some authors (Dole et al., 2011) argued that the summer of 2010 was the result of the internal variability of the climate system, rather than a clear response to global warming.
Many works have analysed the influence of the large-scale dynamics on the variability of extreme temperature indices using several methodologies. For example, Della-Marta et al. (2007) and Carril et al. (2008) related the influence of the atmosphere dynamics and sea surface temperature (SST) to heat waves using methods based on empirical orthogonal functions (EOFs) or canonical correlation analysis (CCA). Other studies use circulation types (CTs; Yiou and Nogaj, 2004;Yiou et al., 2008;Van den Besselaar et al., 2010;Fernández-Montes et al., 2013). These are guided by the results obtained by Corti et al. (1999), who pointed out that recent climate changes can be interpreted in terms of changes in the frequency of occurrence of atmospheric circulation regimes. Some patterns related to anticyclonic and blocking situations favour the development of warm extreme events over Europe (Yiou et al., 2008;Carril et al., 2008;Pfahl, 2014). Therefore, trends in the appearance of these situations could be the cause of the observed trends. Many studies using CTs have provided information in relation to the increase observed in the frequency of such patterns since the second half of the past century (Huth, 2001;Kysely and Huth, 2006;Philipp et al., 2006;Cony et al., 2010;Bermejo and Ancell, 2009;Fernández-Montes et al., 2013;García-Valero et al., 2012). On the other hand, based on the robustness of the evidence from multiple models, the last IPCC report (Stocker et al., 2013) concludes that it is likely that human influence has altered sea level pressure (SLP) patterns globally since 1951. In that way, climate change induces changes in cir-culation that can further modify the occurrence of extreme events.
However, not all studies find strong links between trends in the occurrence of extreme episodes and trends in the frequency of CTs. Jones and Lister (2009), Bermejo and Ancell (2009) and Fernández-Montes et al. (2013), among others, found that trends in extreme temperatures are mainly due to the increase of temperature within the CTs, rather than the increase of the frequency of occurrence of the CTs. This could suggest that trends in extremes can be in addition linked to other forcings, such as global warming, teleconnection phenomena (El Kenawy et al., 2012;Della-Marta et al., 2007), dryness of soil, increase of the SST, etc. Discrepancies in relating both kind of trends (in frequency of the CTs and EHDs) might be due to the way that CT classifications are built. An example of this can be found in Fernández-Montes et al. (2013), where a study on the relationship of trends in extreme T x indices and changes in the frequency of CTs is presented. In this work an established CT classification previously obtained by the authors was used. This classification, as well as the majority of classifications proposed in the literature, is obtained using a huge number of days corresponding to a long period of time (general CT classifications, hereafter). The problem of these classifications for their application to the analysis of extremes is the low statistical load that these specific atmospheric situations, drivers of extreme occurrence, have against the rest of atmospheric situations. This means that these particular situations would be assigned to generalist clusters, since the assignation is controlled by statistical rules that look for a global optimum. Therefore, relating the frequencies of such generalist CTs to extremes might not be appropriate. On the other hand, the use of single composites instead of CTs could solve the above problem, but it has some drawbacks. One is the restriction of the explanation of all extreme occurrences to a unique atmospheric pattern, when there could be several quite different synoptic situations driving to same kind of extreme. Another drawback is the impossibility of its application in attribution exercises since the efficiency would be 100 %, while the experience says that quite similar atmospheric conditions do not necessarily lead to the same effects at the surface, since other forcings could act. Therefore, the use of CT classifications built specifically for the analysis of extremes is the most reasonable methodology for associating large-scale atmospheric conditions to extreme events.
Another aspect to consider is how trends in the frequency of extremes are obtained. Most studies trying to relate local changes (using local series) to changes in dynamics could be problematic for attribution exercises since other forcings, different to dynamics, would have an important control on the variability of these particular events. Regional series composed of a number of local series contain the information of bigger areas (usually of homogeneous orography), filtering out the local noise. Therefore, these are more appropriate for this kind of studies. In principle, series representing a larger area are more homogeneous, and large-scale dynamics has a larger control on its variability. Statistical downscaling is based precisely on this fact, where estimations of local variability are based on predictors, normally large-scale atmospheric fields (representing the dynamics) and the training of statistical models to consider the effects of other forcings which act at more local scales (Wilby and Wigley, 1997). In addition, the larger homogeneity of regional series than local ones is manifested on its extended use in processes of homogenization of climate series. Most methods of homogenization use reference series, formed by the combination of local series to test the relative homogeneity of the local ones (Peterson and Easterling, 1994). In addition, the analysis of frequency of regional series is largely extended on hydrology and water resource applications (Hosking and Wallis, 2005). Therefore, regional series should be used in attribution works in order to reinforce the signal imposed by the dynamics on regional variability. However, it should be checked whether the regional series represent the behaviour of their constituent local ones. Using suitable regionalization procedures should assure this compliance.
This study focuses on mainland Spain and the Balearic Islands for several reasons. First, this region has experienced a significant increase in T x during summer in the last decades Rodríguez-Puebla et al., 2010;Hertig et al., 2010;El Kenawy et al., 2011). Furthermore, trends would have continuity during all of this century considering the results of the AR4 (Stocker et al., 2013), which also provides information about the great sensitivity of the region to future changes, especially in the projections of T x in summer. Second, the complex orography of the study area causes a large spatial variability of climate impacts that the same atmospheric synoptic situation has (Font-Tullot, 2000), making of great interest the studies based on looking for regional differences. In addition, the availability of a highresolution gridded data set developed over the region in the last years (Herrera et al., 2010;Gómez-Navarro et al., 2012) for the analysis of variability facilitates carrying out studies like this.
In this work a novel way of attributing observed trends in the frequency of EHD occurrence to trends in the frequency of CTs is presented. Unlike previous works, here it is proposed a CT classification built based on the definition of extreme (without using general classifications). For the construction of these classifications, regional information contained in the adopted definition of extreme is used, reinforcing on this way the relationships between EHD variability and dynamics. Furthermore, different CT classifications are obtained in order to test the impact that different combinations of atmospheric fields (used for defining the CTs) have on the final results, this being another important difference with respect to other previous works. In addition, the results of the T x regionalization in summer contribute to give new insights supplementing the results of previous works performed over the Iberian Peninsula (IP). This con-tribution is organized as follows: Sect. 2 describes the data sets employed. In Sect. 3 the regionalization procedure necessary for adopting the definition of extreme and the analysis of the obtained regions are presented. Section 4 shows the method followed for characterizing the CTs based on the regional extreme definition as well as a comparative study of the six obtained CT classifications using different combinations of variables. The way of obtaining the links (efficiencies) between the CTs and EHD occurrences, the exercise of attribution and its results, and an analysis of the stability of such links are explained in Sect. 5. Main conclusions and discussions are in Sect. 6.

Surface temperature data
Several high-resolution climate databases for the IP (Herrera et al., 2010) or including it Haylock et al., 2008) have been developed during the last years. These databases have been built by interpolation techniques applied to, in principle, a dense observation network. Although generally reliable, these databases present some known inconsistencies (Gómez-Navarro et al., 2012) due to differences in the raw observational series, the interpolation method, or the different quality controls applied to the data.
Daily T x series of the Spain02 (Herrera et al., 2010) gridded data set are used for this work. This data set was chosen mainly due to the larger number of stations used compared to other similar products available. The large spatial resolution over mainland Spain and the Balearic Islands (0.2 • × 0.2 • ), together with the length of the period , ensures a sufficient spatial and temporal coverage over the study area. Since the work focus on extremely hot days, only dates between 16 June and 15 September are considered.

Large-scale atmospheric data
The data for characterizing the structure of the atmosphere consist of daily fields at 12:00 UTC of SLP, temperature at 850 hPa (T850) and geopotential height at 500 hPa level (Z500) extracted from the ERA40 reanalysis (1958-2002Uppala et al., 2005) and ECMWF analysis (2003)(2004)(2005)(2006)(2007)(2008). The maximum common resolution (1.125 • ) is used for the period 1958-2008. The variables considered are commonly used for the diagnostic of meteorological situations potentially leading to extreme heat events. In this context, SLP offers information about fluxes at low levels and, hence, about the area of provenance of the air mass reaching a given region. T850 provides information about the temperature at low atmospheric levels, tightly related to surface temperature (Brands et al., 2011). Finally, Z500 provides a global vision of the mean atmospheric state. Furthermore, it provides some insight into the overall trough and ridge patterns over the study area, indicating large-scale advection and subsidence in the atmosphere (Sheridan et al., 2012).

Regional series and EHD
This section describes the method of regionalization followed to identify the regions whose constituent local series have similar temporal variability. As was discussed in the introductory section, regional series respond better to changes imposed by dynamics filtering out the signals controlled by other forcings affecting the variability of local series. The definition of EHDs from a regional point of view, the analysis of the variability of the regional series, and the analysis of the coherency of the regions as well as the differences among them are also exposed in this section.

Clustering procedure for regionalization
The procedure applied is similar to that employed by Jiménez et al. (2008) and Lorente-Plazas et al. (2014). First, a Principal Component Analysis (Storch and Zwiers, 1999) in Smode is applied to the correlation matrix calculated using daily anomalies of T x (obtained with respect the seasonal cycle). Only the three first EOFs following the scree plot test (Cattell, 1966) were retained, explaining more than 80 % of total variance. Second a two-step clustering method is applied to the loadings of each grid point in the retained EOFs. The Ward algorithm (Ward, 1963) is employed for obtaining the number of groups and centroids that are used as seeds for a definitive K-means clustering (Hartigan and Wong, 1979). A more detailed explanation of the regionalization method can be found in Lorente-Plazas et al. (2014).

Regions
The aim of this study is not to perform an exhaustive analysis of the regions, but rather to use the regional series as a tool for achieving the attribution goal, so only some aspects related to the time variability of these series are presented below. Figure 1 shows the eight obtained regions. Their names have been established according to their geographical locations: SW, NE, E, Cs (southern central), NWs (south of northwest), NWe (east of northwest), NWw (west of northwest) and N. Regional series have been constructed by averaging the time series of the grid points belonging to the same region. Table 1 (first four columns) shows some statistics of the regional series: mean, trend, standard deviation and 95th percentile. A meridional gradient of mean and percentile values is observed, with the warmest regions being located in the southern half of the country (SW and Cs).
The temporal variability of the regional series is presented in Fig. 2. Solid black lines represent the T x mean seasonal series. For all regions two different periods are observed, confirming the results obtained in previous works that analysed the evolution of summer T x series over the IP (Brunet et al.,  Table 1. Statistics summarizing the eight regional series. The second and third columns show the mean and trend ( • C decade −1 ) of the daily maximum temperature series. The standard deviation of the de-trended series is shown in the fourth column. The fifth column exhibits the 95th percentile. Finally, the sixth column shows the EHD trends (days/decade 2007) and in further Mediterranean regions (Burić et al., 2014). The first lies between 1951 and 1977, when temperatures dropped significantly, with 1977 being the coldest year for most regions. The second period  is characterized by a significant rise of T x . The 1990s were an especially warm decade, with the hottest year occurring for most regions. In particular, 1994 was the hottest for eastern (NE, E and Cs), 1991 for western (SW and NWs) and 1990 for northern regions. During the last decade, a decrease of the standard deviation is observed in the central, southern and Figure 2. Temporal evolution of the eight regional series (see Fig. 1). The seasonal mean series of maximum temperature are represented by the black bold curve (right y axis). The yearly number of regional EHDs and its running mean series (of 11 years) are depicted by vertical bars and the grey curve, respectively (left y axis).
eastern regions, and an increase in the northern ones (N, NWe and NWw).

Extremely hot day definition
The definition of extreme is adopted using the 95th percentile of the regional series (fifth column of Table 1). These percentiles are calculated using the entire period of the Spain02 data set (1951Sect. 2). Thus, a day is defined as an EHD when any of the regional T x series exceeds its 95th percentile. In this way, an ensemble of 863 EHDs were identified in the period 1951-2008. Figure 2 shows the seasonal frequency (bars) and running mean series (11 years) of the regional EHD occurrences (grey line). A noticeable increase is observed in the frequency of EHDs since the 1990s. This finding is in agreement with other studies that found significant increases in climate extreme indices related to T x over the IP Fernández-Montes et al., 2013) and Europe (Klein Tank and Können, 2003). Comparing for all regions the evolution of seasonal T x (black curves) and seasonal frequency of EHDs (bars) of Fig. 2, it is observed that the warmest year does not necessarily coincide with the year of the highest EHD occurrences. In general, the year with most extreme occurrences was 2003, especially in the northern regions (N, NWe, NWs and NE), coinciding with the extraordinary heat wave that affected many countries of western Europe (Beniston, 2004;Trigo et al., 2005). The persistence of EHDs (number of consecutive EHDs) during this year was also extraordinary in most of our regions, reaching up to 10 days in NE, N, NWe and NWs, and up to 16 days in the SW region. In Cs and NWw the maximum number of EHD occurrences took place in the early 1990s, whereas in the E region it curiously occurred in two consecutive years (1958 and 1959).
A trend analysis of the seasonal EHD series has been performed in order to use them for the attribution exercise. Trends and their statistical significance have been obtained by using Sen's algorithm (Sen, 1968) and the Mann-Kendall test, respectively. In addition, to be able to connect these trends to the trends in the frequency of CTs, EHD trends have been calculated for the period 1958-2008 instead of 1951-2008. The election of this period is in accordance with the non-availability of atmospheric data used for obtaining the CTs in the period 1951-1957 (see Sect. 2). The sixth column of Table 1 shows the trends obtained for each region as well as the confidence intervals (95 % of confidence) of these trends. All regions have positive and significant trends (95 % of significance in all regions but 90 % in NWe and E). Largest trends are observed in the inner regions (Cs, NWs and NE), showing a pattern very influenced by the distance to sea, similar to that reported in Bermejo and Ancell (2009) and in Gómez-Navarro et al. (2010), which analysed spatial warming patterns over the IP.
Another important aspect of the obtained regions is its internal coherency and differences among them. Coherency of the regions can be estimated attending to the degree of representation of the region with respect to its constituent series. Thus, a regional EHD should be a signal of many EHD local occurrences at the same time in many of their grids. In this sense, the mean percentage of grid points belonging to the same region that experiences local EHD occurrences when a regional EHD is observed has been obtained for each region. Results show that around 60 % of grid points for most regions have a local EHD when a regional EHD occurs, confirming that regional series, obtained by averaging all the constituent series, are a good thermometer of the global behaviour of the region. NWs and E regions are the most and least homogeneous, with percentages of 65 and 52 %, respectively. In addition, if the 90th local percentile is considered, more than 80 % of grid points overcome this percentile when in the region an EHD occurs. Differences among the regions can be obtained by calculating the probability of simultaneous EHD occurrence among them. Table 2 shows a symmetric matrix with these probabilities. The symmetry of the matrix is in accordance with the number of EHDs in the different regions being the same because of the use of the 95th percentile of each region (all regions have the same number of EHDs during the comparison period, 1951-2008). Diagonal values of the matrix show the probability of occurrence only in this region (without EHD occurrences in the remaining regions). Results indicate that the E region is the most distinctive (40 % of non-simultaneous occurrences), whereas NWs shares many episodes with many of the regions (5 % of occurrences only in this region). The lowest probability of simultaneous occurrence is between E and NWs regions (15 %), and the largest between NWs and Cs (62 %). In general, there are important regional differences, which points to the conclusion that a given large atmospheric pattern could have a different effect over the various regions.

Characterization of EHD circulation types
As discussed above, links between EHD occurrences and CTs could not be well established when general CT classifications are used. We propose the obtainment of specific CT classifications based on the variable to study, in this case EHD occurrence. Before obtaining the final CT classification, main synoptic situations leading to EHD occurrences at the different regions are characterized. In addition, since the results could depend on the atmospheric fields used for characterizing these situations, six different combinations are examined. This will permit testing the sensibility of the final results to the chosen fields.

Clustering procedure for CT characterization
To characterize the atmospheric patterns drivers of EHD occurrences, a similar procedure to those used for general CT classifications is followed. The main difference is that only days characterized as EHDs are used for the clustering. The two-step method applied in García-Valero et al. (2012) is used here. In the first step, a principal component analysis in T mode (PC-ModeT) clustering is performed (Kysely and Huth, 2006). This clustering allows obtaining the necessary seeds for the second clustering step, defining the number of clusters (a priori unknown). Second, a K-means algorithm is applied over the retained PCs, using for initializing the clustering the seeds obtained in the previous step. The geographical window for clustering is identical to that employed in García-Valero et al. (2012), which covers completely the IP and Balearic Islands (35-45 • N and 10 • W-6 • E). Performing the clustering using larger windows might lead to inclusion of more noise in the classifications, obtaining probably some clusters representing further dynamical structures with little influence on the regional climate variability of the IP. Despite the windows employed covering a small area, the results of clustering are representative of bigger areas (García-Bustamante et al., 2012), so a larger window is used for representation of the centroids (composites of the clusters), allowing for better visualization of the synoptical structures. Six CT classifications have been obtained. Three of them consider the atmospheric fields individually (SLP, T850, Z500). The other three consist of the combination of all possible pairs of fields (SLP-T850, SLP-Z500 and Z500-T850). The number of clusters of each CT classification depends on the retained PCs used for the PC-ModeT clustering, this being double the retained PCs (García-Valero et al., 2012). For classifications with only one atmospheric field, six clusters are obtained, whereas eight clusters are obtained if two fields are considered. Therefore, three or four PCs are retained if fields are considered individually or in pairs, respectively. In all cases the explained variance by these PCs is over 90 %.

Evaluation of the CT classifications for EHD description
The way of selecting those days defined as EHDs (any EHD occurrence over almost one region) makes an ensemble of situations affecting the various regions in a different way. For a given region, some of the CTs will have great influence    on its EHD occurrences, but some others do not have (or little) influence, affecting other regions more. On the other hand, using clustering techniques for classifications, as well as to force the clustering of a great number of events to a reduced number of clusters (6 or 8), causes the existence of noise in the classifications. A question arising is if there is any classification, among the six obtained here, that characterizes better the EHD occurrence for each of the regions. In order to address this question, an index, effectiveness index (EI), is defined. This index is calculated as the ratio between the number of EHD occurrences and non-occurrences in the region under a given CT. Therefore, for a given CT classification and region, there is a set of EIs. The best classification is chosen from that whose standard deviation and range of the EIs are the largest, since it separates better the most-from the least-influential CTs. In addition, those CTs with EI > 1 for a given region are defined as "extreme CTs", because the probability of extreme occurrence is larger than the non occurrence. Table 3 shows the range and standard deviation (range-SD) obtained for all the regions and CT classifications. Last column shows for each region its best (highlighted in black) CT classification with its extreme CTs. Results indicate that two synoptic variables characterize better the EHDs for most regions; therefore, hereinafter the remaining analysis will focus only on the CT classifications formed by two atmo-spheric fields. Z500-T850 is the best one for NE, NWw and Cs (inner-eastern regions and northwestern). SLP-Z500 characterizes better the E region. SLP-T850 is the best for SW, NWe and N (western and northern regions). T850 classification is the best for the NWs region. However, in this area SLP-T850 performs quite similar. Hence, and for the shake of clarity in the analysis, SLP-T850 has been considered as the best classification in this case.
Figures 3-5 show the centroids (composites) of the CTs (contours) belonging to the CT classifications formed by two atmospheric fields as well as their efficiencies (shaded) over the different regions. Efficiency is defined as the conditional probability of having an EHD in a region under a given CT. Another interesting parameter is the contribution of a given CT to the occurrence of the EHDs in a given region. The contribution is assessed by calculating the ratio between the number of observed EHDs under a CT and the total EHDs observed in the region. Tables 4 and 5 depict the efficiencies and contributions values.
The efficiency patterns are quite similar for the three CT classifications. However, for each region the efficiency is higher for the classification that gives larger spreads in the EI (Table 3). Some examples follow. The efficiency pattern related to CT8 is equivalent in all classifications and shows high efficiency over the E region. The efficiency is higher for the SLP-Z500 classification, which has the largest EI for E Figure 3. Centroids of the CTs for the SLP-T850 classification. SLP and T850 are represented by contours, black lines for SLP and colour lines for T850. Shading denotes the regional EHD efficiencies associated with each CT. The top left corner shows the number of the represented CT. Right corner indicates the frequency (in percentages) of each CT, before and after applying the allocating method (Sect. 5.1).
region. Similar results are found for CT5 of SLP-T850, CT6 of Z500-T850 (best) and CT5 of SLP-Z500 in the NWw region, and for CT3 of SLP-T850 (best), CT3 of Z500-T850 and CT1 of SLP-Z500 over the SW region. These results suggest that some CTs belonging to different classifications are equivalent; i.e. they give similar efficiency patterns. In fact, this can be corroborated by calculating the common days of the mentioned CTs (not shown). This highlights the need for studies on the sensitivity of the CT classification to the atmospheric variables employed.
The comparison of the atmospheric situations, among the different classifications, associated with similar efficiency patterns, enables some conclusions to be drawn about the main drivers related to the EHD occurrences. Regarding the T850 variable, the regional efficiency shows the highest values in regions where temperatures are near and above 20 • C. This feature is very common in many CTs of the different classifications such as the CTs 2, 3, 4, 6 and 8 of SLP-T850 and the CTs 1, 2, 3, 4, 5, 7 and 8 of the Z500-T850 classification. The wind provenance, inferred considering the SLP field, is also an important factor contributing the occurrence of EHDs in some regions, because of the warm advection over specific regions. Many CTs are related to situations of wind blowing from inner towards coastal areas, causing the highest efficiencies in the latter. The inner IP regions are highland plateau areas where high temperatures are observed (see Table 1). When wind blows from this area towards the sea through valleys, air is adiabatically compressed, causing an important warming in lowland regions. Some examples of these situations can be identified in the classifications by analysing the efficiencies and contributions of some CTs. Five regions are mainly affected by this: NWw, especially when wind blows from the east because of the presence of high pressures over western Europe (CT5 of SLP-T850 and SLP-Z500); SW in northeastern wind conditions (CTs 1 and 3 of SLP-T850 and CT1 of SLP-Z500) as a result of the presence of high pressure over the Mediterranean and relative low pressures over the southwest of the IP; the E region under strong western zonal wind (CT8 of SLP-T850 and SLP-Z500), induced by the location of high and low pressures over the Atlantic and Mediterranean, respectively; and the NE and N regions in southwestern wind situations (CT6 of SLP-T850 and SLP-Z500). Conversely, CTs with weak SLP gradients (stagnant situations linked to thermal lows) are specially important in the EHD occurrences at Cs and NE (CT4 of SLP-T850 and SLP-Z500), and N, NWe and NWs (CT2 of SLP-T850) regions. Such situations were also pointed out by Pfahl (2014) like those more relevant guiding hot extremes in summer over the IP. Regarding the Z500 field, EHDs are associated with large amplitude ridges over the IP. Regions with the highest efficiencies are located to the west of the ridge axis, where there are high stability and warm advection. The efficiency and contribution are only important in the E region when zonal wind at 500 hPa is strong (CTs 5 and 8 of Z500-T850 and CT8 of SLP-Z500). This pattern usually takes place when hot episodes over the IP are ending up. The analysis of the CT transitions (not shown) for consecutive EHD episodes reinforces this result. The hottest areas travel from western to eastern regions, following the ridge movement.

Linking EHD trends to CTs
Significant and positive EHD regional trends were found for all regions (Sect. 3.3). An interesting question is whether there is any relation between these trends and changes in the frequency of occurrence of the CTs. Therefore, the existence of trends in the frequency of CTs should be assessed. The main problem is that frequencies of CTs only can be ob-tained during extreme episodes. There are atmospheric situations, similar to that described in the centroids, that were not included in the previous clustering step because EHDs did not occur in any region. As was discussed above (Sect. 1), the same atmospheric patterns could lead to some different weather situations. To take this into consideration, the rest of the days initially not considered for characterization have to be assigned to one of the obtained clusters. For this reason, links (efficiencies) between CTs and EHD occurrences have to be recalculated. The method of assignation of such days is presented below.

Allocation method
Each CT is the result of a set of similar atmospheric conditions being its associated centroid the mean value (composite) of this ensemble or population. However, there is some dispersion inside each group. To allocate an atmospheric situation (not considered as an EHD) into the CTs means finding the population where such situation fits better. For this task the election of some metrics is fundamental. Two metrics have been used for assignation: the Spearman correlation and the Euclidean distance. The former measures the spatial similarity of two fields (Hofer et al., 2012;Vautard Figure 5. Same as Fig. 3 but for the SLP-Z500 classification. SLP and Z500 are represented by isobars (black lines) and contour levels (blue shaded lines). and Yiou, 2009), whereas the latter evaluates the differences in the intensity of patterns (Fettweis et al., 2011). Both metrics are relevant and complementary (Nogaj et al., 2007;Fettweis et al., 2011). The population of a given CT is characterized by the distribution of the distances to and correlations with the centroid using only the days considered in the characterization (Sect. 4). To obtain the population of each CT, the atmospheric data are first standarized by grid point; then correlations and distances to the centroid are obtained. Mean and standard deviation for standardization are obtained using all days (extremes and non-extremes) of the study period. A criterion to allocate a situation into a given CT is that the distance to (correlation with) the centroid be lower (higher) than some given thresholds. Hence, the election of these thresholds is a problem to solve. On the one hand, thresholds could be fixed to certain percentiles of the populations, but this could be problematic because this would have strong subjectivity. On the other hand, thresholds could be the maximum (minimum) values of distance (correlation) obtained from days clustered for the characterization of the CT leaders to EHD occurrences (Sect. 4). We have chosen the latter option because this is free of subjectivity. shows the thresholds considered for all CTs and classifications. Some thresholds of certain CTs are too broad, allowing the assignation of more days into them. These CTs are those with higher noise inside them, and probably they are formed by a larger heterogeneity of situations. By contrast, there are CTs more homogeneous with thresholds more restrictive. In general the most homogeneous CTs agree with atmospheric patterns having high efficiency in most regions simultaneously. Finally, a given day is allocated to the CT on which distance and correlation are lower and higher, respectively, than the thresholds defined for it. Days that can not be assigned to any CT are then allocated to a new group (unclassified group, CT9; Seubert et al., 2014). This way of assignation might derive from the allocation of a situation to more than one of the established centroids. In these cases, it is allocated to the nearest cluster considering only the Euclidean distance. The use of this criterion is in accordance with the lower spread of the Euclidean distance populations than that observed for correlation. Table 7 shows the number of days belonging to the different CTs before and after the allocation process. Approximately 70 % of the days for all classifications are assigned, whereas the rest are allocated to the unclassified cluster (CT9), with Z500-T850 (SLP-Z500) being the one with most (fewest) unclassified days (33 % vs. 28 %). There is a large variability in the increase of the days belonging to the different CTs. CT1, CT6 and CT8 have the largest increase (10 times) in all classifications, whereas others like CTs 2 and 4 of SLP-T850 and CTs 2 and 4 of Z500-T850 present small changes (less than 20 % of the initial clustered days). The assessment of the quality of the clusters before and after the allocation is of major relevance. In order to ensure the reliance of the assignation method, the explained cluster variance (ECV) of each classification (Table. 8) is analysed. Results show few dif-ferences among the classifications. Z500-T850 classification has the best quality, observing even an increase in the quality of clusters after the allocation. The quality of the other classifications worsens slightly after the allocation, with the SLP-T850 classification being slightly better than the SLP-Z500 one. On the other hand, when comparing the quality for the three classifications after the allocation (ECV > 46 % in all cases; see Table. 8) with some others obtained using general CT classifications, even better results are observed in our classifications. Hence, in Philipp et al. (2006) an ECV lower than 40 % was obtained in summer, and in García-Valero et al. (2012) of 47.2 %. These results support the suitability of the method followed for allocation.

Analysis of the allocation
It is instructive to explore the consequences after applying the allocation procedure on the properties of the groups or CTs. One consequence of the larger number of classified days is the decrease in the efficiency of the CTs (see Table 4). Obviously, the lower the increase of the number of days, the smaller the decrease of the efficiency. This effect stands out in CTs 2 and 4 for the SLP-T850 and Z500-T850 classifications, which have the highest efficiencies in most regions. Changes in the efficiency after the allocation also affect the EI used for deciding the best CT classification for each region (Sect. 4.2). Following the same criterion as above, results remain unaltered for most regions, except for E and NWw regions, which are now better characterized by the SLP-T850 classification.
The shape of the populations of distances and correlations of each CT can be also affected. The net effect of the allocation process is to include new days located further from the centroid. Nevertheless, a small number of clusters hardly change the populations after the assignation (CTs 2 and 4 of SLP-T850, CTs 2 and 4 of Z500-T850 and CTs 4 and 7 of SLP-Z500), which are coincident with the ones with higher efficiency. Figure 6 shows two examples of the correlation histograms before and after the allocation process (top) as well as the empirical cumulative distribution function (bot-   q q qq q q q q q q q q q q q qq q q q q q q q qq q qq q q q q q q q q q qq q q qq q q q q q q q q q q q q q qq q qq q q q q qq q q qq q q q qq q q q q qq q q q q q q q q q q before after q q qq q q q qq q q q q qq q q qq q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q qq q qq q q q q q q q q q q q q q qqq q q q q q q q q qqq q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q qq q q q q q q q q q q q qq q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q qq q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q tom). The left and right panels show an example of great and small changes in population, respectively. Once the allocation is performed, an analysis of trends in the frequency of CTs is carried out. Table 9 shows the trends obtained for each CT of the different classifications. Trends and their statistical significance are obtained following the same methods (Sen's algorithm and Mann-Kendall test for trend and its statistical significance, respectively) as used for estimating the regional EHD trends. Results indicate that only two CTs of each classification have trends with statistical significance (one more but without significance). Considering this and their associated efficiencies (Table 4), it could be said that trends in CTs would be linked to the increase of EHD occurrences in northern regions (mainly N and NE). It is interesting also to observe that the CT9 (formed by the unclassified days) has a negative trend for the three classifications, meaning that in time a larger number of days are assigned to the centroids linked to EHD occurrences.

Attribution of EHD trends
Temperature changes can be linked to several factors, one of them being the changes in atmospheric circulation. Trends in the EHDs could be considered also as an indicator of temperature changes. In this subsection, a model to attribute regional EHD trends to trends in the frequency of CTs is presented.
The regional EHD trend in a given region, T r , can be written as the sum of two terms: where T r c is the trend attributable to changes in atmospheric circulation and T r o the trend related to other factors. Now, T r c can be described as a linear function of the changes in the frequency of the CTs. We propose the simple model where the subscript i denotes the CT number, from 1 to n (the number of CTs); T i is the frequency trend of the CT i ( Table. 9); and r i denotes the efficiency of the CT i over the region r (Table. 4), calculated after the allocation process.
Using this simple model T r c can be calculated for all regions using the different CT classifications. The results are summarized in Table 10. The last column depicts the observed regional EHD trend. With independency of the CT classification, the reproduced trends are in all cases smaller than those observed for a given region, but they are strongly dependent on the CT classification, except for SW and NWw regions where few differences are appreciated among classifications. In general terms, the Z500-T850 classification (that of the best quality; see Sect. 5.2) obtains the highest reproduced trends except for NWw. Using this classification a fraction between 40 and 50 % of the observed trends is attributed to changes in the frequency of CTs in the E, NE, NWe and N regions, about 20 % in the western and southern regions (SW, Cs and NWs), and only 14 % in NWw.

Within-type variations in the efficiency
In the exercise of attribution, links between CTs and EHD occurrences, established by the efficiency, have been taken as long-term mean values of the efficiency throughout the entire study period , because our main goal is to attribute long-term regional EHD trends. On the other hand, if our purpose were to describe the variability at higher frequencies, an analysis of the stability of links (obtaining the links using higher frequencies) should be carried out. The instability of links is frequently known in the literature as within-type variations (WT). Some of the limitations of the use of CTs for downscaling purposes is precisely the existence of WT, meaning that other factors could be controlling the links. For example, the scarcity of precipitation in spring would lead to drier soils, and during summer landatmosphere processes would be more intense (Seneviratne et al., 2010;Jerez et al., 2012), enhancing air temperature and efficiency. Another example is the persistence of atmospheric situation drivers of high T x ; this also would cause an enhancement of land-atmosphere processes. Persistence might be originated by the existence of SST anomalies caused by teleconnection phenomena . In order to evaluate the possibility of using the CTs obtained here for downscaling of annual EHD in the different regions, an analysis of the stability in the efficiency for the different regions has been carried out. To do so, annual efficiency series have been obtained, and, in a way similar to other works (Jacobeit et al., 2009;Fernández-Montes et al., 2013), running mean series of 31 years are derived for the analysis of the stability. Figure 7 shows the obtained running mean series for every region and for the CTs composing the best CT classification for each region (that with the highest dispersion of the EI after the allocation process; see Sect. 5.2). In general terms, those CTs with the highest efficiencies in each region increase its efficiency by 10 % throughout the study period . This is the case of CT2 in the regions Cs, NWs, NWe, NWw and N, and CT4 in the E region. Anyway, there are some other important WT variations in CTs affecting the inner regions more, pointing to a rise of the efficiency. Some examples follow: CT-3/4 for Cs, CT-1/2 for SW and CT-1/3/7 for NWs. Another interesting result is the increase of the efficiency of CT4 in the northern regions (NWw, NWe and N), while it decreases in the E region, confirming once again the remarkable regional differences existing over the studied area. Among all the regions, the least affected by WT is the NE region, where only CT1, with low efficiency over the region, increases its efficiency. This last result is curious because this region was noted by Bermejo and Ancell (2009) as one of the Spanish regions where larger WT in T x were observed during the second half of the past century, observing warming in T x inside many of the CTs obtained by these authors. Differences between our results and those found in Bermejo and Ancell (2009) could be in the variable analysed, in our case extreme events and in their case T x mean. Hence, the increase of T x mean may not mean an increase in the efficiency of EHD occurrence.
Some reasons for the increase in the efficiency found in the CTs can be linked to the increase of the persistence of CTs. To evaluate this possible influence, the correlation between the de-trended seasonal frequency series and mean seasonal persistence series of the CTs with the most important WT have been calculated. Results show positive and significant correlations, between 0.6 and 0.7, which support the influence of the persistence on the efficiency rise. In addition, the decline of soil moisture observed over the IP since the 1970s (Sousa et al., 2011) could be another factor contributing to higher efficiencies.

Conclusions and discussions
This work sets as its intention the attribution of regional trends in EHD occurrences observed at eight Spanish regions (during the period 1958-2008) to changes in the frequency of atmospheric situations linked to the occurrence of such events. The study is centred in summer, when EHDs have a larger relevance. The method followed is based on the procurance of CT classifications using regional information contained in the same definition of extreme. This differs from others studies based on general CT classifications. Therefore, a regionalization of daily T x series is carried out, deriving from this exercise the regional series used for adopting the definition of extreme. In addition, the internal coherency of the regions and differences among them are also analysed. Later, a characterization of the atmospheric situations, drivers of EHD occurrences in the different regions, has been obtained, centring these classifications only on those days defined as EHDs. In this way, six CT classifications each using a different combination of atmospheric fields have been obtained, and an analysis for finding the most suitable classification for each region has been performed. Links between CTs and EHD occurrences have been defined by means of the efficiency that a given CT leads to a regional EHD in each region. To establish such links, the CT classifications obtained in the characterization step have been extended to the rest of days (non-EHDs). To do so, a method of allocation of these days to the centroids obtained in the characterization is presented. Finally, the method of attribution is presented as well as an analysis of the stability of the links between large-scale dynamics and EHD occurrence. The main conclusions are summarized below.
-Eight regions with different daily T x variability are identified. Regions have high internal coherency and important differences among them. Positive and significant regional EHD trends are found across most regions. Generally, such trends are larger in central and northern regions, and lower in the SE and NWe regions.
-To find CT classifications describing the extreme occurrences, it is convenient to base them on the variables to analyse, as well as on a characterization of the associated atmospheric patterns using exclusively those days defined as extremes. The method proposed here produces CT classifications of similar quality to those obtained in general CT classifications which use other clustering techniques. However, this method ensures a more precise allocation of the extreme days to the correct clusters than general CT classifications.
-The choice of the most suitable combination of atmospheric variables to define the CTs becomes of major relevance. This study finds that classifications using combinations of two atmospheric variables generally perform better than those using only one. In the case of the studied area, SLP-T850 characterizes better the EHD occurrences for most regions. However, other combinations, such as Z500-T850, have better results in the Cs and NE regions. In addition, Z500-T850 is the CT classification with the best quality, whereas SLP-Z500 has the worst quality in terms of the explained variance by the clusters.
-In the different obtained classifications, only a small number of CTs have significant trends of its seasonal frequency, a result which is in line with other studies using general CT classifications (Fernández-Montes et al., 2013;Bermejo and Ancell, 2009). Such trends depend on the atmospheric variables used for the classification, with the largest trends occurring in the classification of the best quality (Z500-T850). CT2 of the Z500-T850 and SLP-T850 classifications have the largest trends. This atmospheric pattern is associated with high occurrence of EHDs in most regions simultaneously, specially in the central and northern regions.
-Part of the EHD trends observed in some regions can be attributed to changes in the CT frequencies. However, the trends reproduced by the attribution method are in all cases lower than the observational ones, indicating that part of the trend has to be attributed to other factors. The attributed trends have a great dependence on the region, as well as on the CT classification. Thus, the bestquality classification, Z500-T850, is able to attribute the larger part of the observed regional EHD trends for most regions except for NWs. Considering this classification, a fraction between 30 and 50 % of the observed trends is attributed in the E, NE, NWe and N regions, about 20 % in the western and southern regions (SW, Cs and NWs), and hardly 14 % in the NWw.
-Some of the most-influential CTs in terms of leading to EHD occurrences in the different regions present WT variations in relation to their efficiency, except for NE. Approximately the most efficient CTs have increased their efficiencies by 10 % during the study period. These changes may be associated, among other things, with the increase of the persistence of such CTs through time.
The attribution exercise reveals that the observed regional EHD trends can be only partially attributed to changes in the atmospheric dynamics and that they have an important regional component. This suggest that there are other factors involved in the EHD trends -such as global warming, soilatmosphere feedbacks or changes of surface properties -that contribute to increasing positive changes in EHD frequency. A fact that reinforces this asseveration is that the trends in central regions are less affected by changes in atmospheric dynamics. This could be related to a depletion of soil moisture that enhances the positive land-atmospheric feedback, as previously stated by other authors . This also supports the WT variations mentioned above.
The links obtained here between regional EHD and CTs, and the allocation method followed, can be used for estimating the role of the atmospheric dynamics in the long-term regional changes of EHD occurrence over the different regions under different climate change scenarios. This is a natural extension of this work. In addition, the methodology applied here could be extended to other extreme events such as floods, droughts, heat waves, etc. The robustness of the regional series as well as the identification of the best largescale atmospheric variables characterizing such events could be of crucial importance when trying to relate regional extreme behaviour to atmospheric dynamics.