1 Introduction

Passive microwave satellite data reveal a record surface snowmelt extent over the Greenland ice sheet (GrIS) during the 2007 and 2008 summers (Tedesco et al. 2008a, b), in line with the decadal melt trend that has been observed and simulated since 1990 (Mote 2007; Fettweis 2007; Hanna et al. 2008). At the same time, Arctic sea ice cover reached an absolute minimum extent in September 2007 (Serreze et al. 2007; Comiso et al. 2008), while GrIS iceberg discharge (Howat et al. 2007) and glacier velocity have also increased (Rignot and Kanagaratnam 2006; Rignot et al. 2008). While these observations may not reflect long-term variations, because of the considerable year to year variations in GrIS mass balance (Howat et al. 2008; Wouters et al. 2009; Van den Broeke et al. 2009), all studies show a clear trend towards GrIS mass loss over the last decades (Hall et al. 2008; Thomas et al. 2008). This trend towards enhanced melting is a direct consequence of Arctic warming which has been attributed to increased atmospheric greenhouse gas concentration (Fettweis 2007; Stroeve et al. 2007; Hanna et al. 2008, 2009;). However, abnormal atmospheric circulations (bringing warm air masses to the GrIS) are thought to explain part of the extreme surface melt events of the last years (Hanna et al. 2009; Stroeve et al. 2008; Rogers et al. 2004; Tedesco et al. 2008a). Atmospheric circulation itself is likely to be affected by anthropogenic forcing, resulting in changes in large-scale circulation states.

Mote (1998a, b) was first to link daily melt extent (derived from microwave satellite data) to mid-tropospheric circulation variability during the period 1979–1989. We extend here the study to the period 1958–2009, using output of the regional climate model MAR (Modèle Atmosphérique Régional). The MAR model (Gallée and Schayes 1994) previously demonstrated its ability to realistically simulate surface melt over the GrIS (Fettweis et al. 2005, 2006; Tedesco et al. 2008a); MAR can directly calculates meltwater run-off rate (Lefebre et al. 2003) rather than just melt extent. The simulation of Fettweis et al. (2008) was extended to the period 1958–2009, using the same model set-up as (Fettweis 2007) with some improvements in the radiative and snow model. These improvements reduce the simulated positive bias in melt reported by Fettweis et al. (2008). In order to study the correlation between surface melt and the general circulation over the GrIS, we developed an automatic typeification of daily 500 hPa geopotential height surfaces, based on reanalyses (1958–2009). After describing the reanalyses data in Sects. 2, 3 described and evaluates the Circulation type classification (CTC). The link with surface melt anomalies as simulated by MAR is then discussed in Sect. 4. For the classification, we do not use 500 hPa geopotential height surfaces simulated by MAR because they are similar to those simulated by the reanalyses (used as boundary conditions in MAR).

2 Data

The 500 hPa geopotential height data used here are from:

  • the ERA-40 reanalysis (1977–2002) (Uppala et al. 2005) and the ERA-INTERIM reanalysis (2002–2009) from the European Centre for Medium Range Weather Forecasts (ECMWF). The ECMWF reanalysis is also used to initialize meteorological fields at the beginning of the MAR simulation and to force the lateral boundaries every 6 h during the simulation. The ERA-40 reanalysis is available from September 1957 onwards with a time resolution of 6 h and a spatial resolution of 1.125 latitudinal degree.

  • the Reanalysis 1 from the National Center for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) (Kalnay et al. 1996). The NCEP–NCAR reanalysis is available from 1948 with a time resolution of 6 h and a spatial resolution of 2.5 latitudinal degree.

The ECMWF reanalysis has the advantage of being available at a higher spatial resolution than the NCEP/NCAR reanalysis, and has proven its value over the GrIS (Hanna et al. 2005, 2008). However, changes to the ECMWF model configuration/resolution used for ERA-INTERIM (2002–2009) will negatively affect the homogeneity of the time series. That is why we use both ECMWF and NCEP/NCAR reanalyses to study the mid-tropospheric variability.

According to Fettweis et al. (2008), the yearly GrIS meltwater run-off rate mainly depends on conditions from 1 June to 31 August (JJA). By conditioning surface albedo (Mote 2003), winter accumulation also plays a role in the interannual run-off variability but to a lesser extent than the JJA near-surface temperature. That is why we have limited our CTC to JJA (92 days) over the common period covered by both reanalyses. This yields 4,784 summer days to be classified over 1958–2009. As daily value, we take the 12 h UT values of the 500 hPa geopotential height on a latitude–longitude grid. To give the same weight to all grid points, the data are projected on a 100 km-regular grid. This grid is centred on 72°N, 40°W by using a Lambert Conformal Projection. The grid dimension is 1,400 km by 2,700 km.

3 Circulation type classification

3.1 Background

Several manual or automated CTC approaches have been developed in synoptic climatology; each of them has advantages and drawbacks as shown by the ongoing European project COST733 (Philipp et al. 2010). The manual typeification (often constructed for European or US regions) are highly subjective (Yarnal et al. 2002) and cannot directly be applied to Greenland at the difference of the automatic approaches. The most common automated techniques are correlation-based analyses (Lund 1963; Kirchhofer 1973) and eigenvector-based analyses (principal components analysis (PCA) and empirical orthogonal functions (EOFs)) according to Yarnal et al. (2002). Both techniques have been used successfully to study precipitation and temperature variabilities (Yarnal et al. 2002). Mote (1998b) performed a cluster analysis on component scores of a PCA of daily 700 hPa heights to study the relationship between surface melt and atmospheric circulation over the GrIS.

Here, we propose to use a correlation-based approach to study GrIS melt. In a correlation-based map-pattern classification, each pair of days is first characterized by a similarity index varying generally between −1 and 1. This index is then used to group the days with similar circulations (Lund 1963).

The main differences between the present study and the ones of Lund (1963) and Kirchhofer (1973) are the definition of the similarity index as well as improvements in the final step (the classes construction) aimed at reducing the within-type variability. In addition, our CTC is based on the 500 hPa geopotential height surfaces (called Z 500 hereafter) rather than on the sea level pressure surfaces. This was done because

  • the topography of the GrIS reaches up to 3,250 m at Summit;

  • we show in the next section that the mere presence of an anticyclone at 500 hPa (suggesting a warmer atmosphere) is more important than the flow direction for studying surface melt. We therefore modified the similarity index to take into account the height of the Z 500 surface rather than its pattern, as done in the Lund and Kirchhofer classification schemes.

3.2 2m-temperature versus mid-tropospheric correlation

The melt rate is closely related to the near-surface air temperature according to Ohmura (2001). Figure 1 shows the correlation between the daily mean GrIS 2m-temperature [i.e. a proxy of the surface melt rate as used by Fettweis et al. (2008)] and the mid-tropospheric state. Figure 1a emphasizes that, in both reanalyses, warm events are associated with the presence of an anticyclone at 500 hPa centred above the ice sheet. Such an anticyclone:

  • implies a larger thickness in the lower troposphere suggesting a warmer atmosphere;

  • induces a clockwise circulation above the ice sheet with zonal (west-east) winds over the north of the ice sheet (see Fig. 1b) and southern winds along the west coast (see Fig. 1c).

  • induces downward motion (subsidence) and heating.

The correlation coefficients are however much larger in Fig. 1a than in Fig. 1b and c which reflect the flow directions induced by an anticyclone at the centre of the ice sheet. This suggests that the mean GrIS 2m-temperature depends more on the height of the 500 hPa-(anti)cyclone at the centre of the ice sheet than on its precise location and pattern. Thus, contrary to Lund (1963), the Z 500 surface height takes a larger part in our CTC (via the similarity index) than its pattern.

Fig. 1
figure 1

Correlation between daily mean GrIS 2m-temperature and a the 500 hPa geopotential height, b the 500 hPa meridional wind and c the 500 hPa zonal wind . The ice sheet mask used to compute the mean GrIS 2m-temperature is shown by the black line. The 12 h TU outputs from respectively the ECMWF reanalysis (top) and the NCEP/NCAR reanalysis (below) are used as daily values in the time series

3.3 Algorithm

Let Z 500(d) be the 500 hPa geopotential height surface at 12 h TU at day d.

  1. 1.

    To give the same weight to all grid points in the similarity index computation, the Z 500 surfaces are first normalized, so that the standard deviation is equal to 1. For each grid point (i, j), we have then

    $$ \forall i \in [1,N], \forall j \in [1,M]: \sqrt{{\frac{1} {n-1}}\sum_{d=1}^{n} {(Z_{500}(i,j,d)-\overline{Z_{500}}(i,j))}^2} = 1 $$
    (1)

    where N and M are the grid dimensions, \(\overline{Z_{500}}\) is the temporal mean 500 hPa geopotential height and n (=4,784) is the number of JJA days from 1958 to 2009.

  2. 2.

    For each pair of days, a similarity index based on the absolute difference between the two Z 500 surfaces is computed:

    $$ i(\text{day}_1,\text{day}_2) = 1 - {\frac{1}{c}} . {\frac {\sum_{i,j=1}^{N,M} | Z_{500}(i,j,\text{day}_1) - Z_{500}(i,j,\text{day}_2) |}{N \times M}} $$
    (2)

    where we have chosen c = 2 (The choice of the value of c is discussed in the next section). This index i reaches 1 if the Z 500 surface is compared with itself and is negative if both Z 500 surfaces are very different. Consequently, the closer this index is to 1, the more similar the Z 500 surfaces are.

  3. 3.

    The Lund (1963)-like construction of classes works as following:

    • Each pair of Z 500 surfaces is characterised by a similarity index i ≤ 1.

    • For a given similarity index threshold i c (e.g. i c  = 0.7), we compute, for each day day1, the number of similar Z 500(day2) surfaces, i.e. the number of days day2 for which \(i(\text{day}_1,\text{day}_2)\,{\geqslant}\,i_c.\)

    • The day with the highest number of similar surfaces is the reference surface/day for the first class. This first class contains all days that are similar to this reference day.

    • The classified days (in the first class) are removed from the sample. With the unclassified days, we construct the second class in the same way and so on until all the days are classified.

      In short, the class C k (i.e. a synoptic circulation type) contains all days d such that:

      $$ i(\text{day}_k,d)\,\geqslant\,i_c $$
      (3)

      where day k is the reference day of class C k .

      If we want to classify all days of the sample,

    • either the threshold i c is chosen by the user and the number of used classes increases as i c is near 1. If i c equals 1, the number of classes will be the number of days in the sample. If i c is small (i.e. near zero), a few classes only will be sufficient to classify all days.

    • or the number of classes K is fixed by the user (as is the case here) and the threshold i c is progressively decremented by the algorithm from 1 with a value of 0.01 until all sample days are classified. The smaller the number K, the smaller the threshold i c used to classify all days, and the less similar the days within a single class will be.

  4. 4.

    If the number of classes K is fixed by the user, it is better to use a threshold i c decreasing with the class rank rather than a fixed threshold, in order to avoid most of the days being classified in the first classes (as seen in the next section). This decreasing threshold is defined as:

    $$ i'_c(k) = i_c - \epsilon . ( k - 1 ),\ k \in [1,K] $$
    (4)

    where k is the rank of the class C k and ε is a factor discussed in the following item. This means that the similarity index threshold used to construct the second class is lower than the one used to create the first class and so on until the last class (k = K), where the threshold i c (k) often reaches ∼0.5.

  5. 5.

    Following Buishand and Brandsma (1997), the mean-squared-error skill score S (percentage of explained variance) is defined as as:

    $$ S = 1 - {\frac{MSE_{\text{class}}}{MSE_{\text{tot}}}} $$
    (5)

    where

    $$ MSE_{\text{tot}} = \sqrt{{\frac{\sum_{i,j=1}^{N,M} \sum_{d=1}^{n} (Z_{500}(i,j,d)-\overline{Z_{500}(i,j)})^2}{N \times M \times n}}} $$
    (6)
    $$ MSE_{\text{class}} = \sqrt{{\frac{\sum_{i,j=1}^{N,M} \sum_{d=1}^{n} (Z_{500}(i,j,d)-\overline{Z_{500}k(d)(i,j)})^2}{N \times M \times n}}} $$
    (7)

    and \(\overline{Z_{500}k(d)}\) is the temporal mean Z 500 surface over the days d belonging to the class C k . The CTC quality is linked to the degree of similarity into each class. Therefore the closer S is to 1, the smaller is the variance in the classes (compared with the total variance) and the more similar the days in the classes. The skill score S is very useful to compare the CTCs as well as to assess the 'best’ parameters values to be used in the CTC.

    The value of ε used in i c is computed by the algorithm in order to minimise the skill score S. The algorithm first makes a classification for ε = 0 and increase ε with increments of 0.005 until a absolute minimum of S is reached. This last step takes a lot of computer time since, for each "test” value of ε, the algorithm has to decrease i c until all the days are classified. This last step (made for reducing within-type variability) is original compared to the method of Lund (1963).

3.4 Results

Figure 2a shows the results of our CTC for a fixed number of classes (8) and on a domain centred on the GrIS, using the ECMWF reanalysis data set. The choice of the class number and domain is discussed in the next section. The mean GrIS summer 500 hPa-circulation is a west to east circulation with a south-west to north-east component above the GrIS (see Fig. 4).

Fig. 2
figure 2

The Z 500 average over the days categorized into the eight classes minus the Z 500 average over the sample 4,784 days (shown in Fig. 4). For each circulation type, the proportion of days contained in a class compared to the 4,784 days of the 1958–2009 summers is indicated in brackets. The CTC was built using a the ECMWF analysis and b the NCEP/NCAR reanalysis. Units are metres

The days with a circulation similar to the climatology (Fig. 4) are classified in Class 1 and account for 28% of the total. Type 1 has its highest frequency in July (see Table 1). Types 2 and 5 show negative anomalies above the ice sheet and are associated with a southward displacement of the low normally located over Northwest Greenland. These classes (accounting for ∼33%) are more prevalent in June and in August. Types 3 and 7 (accounting for 24%) display an opposite situation, with a Azores Anticyclone farther north than normal, favouring warm air advection. July, the warmest month of the year, is dominated by these events in addition to Type 1 events. The first three types are the most frequent, but they are also the closest to the climatology. Type 5 (respectively Type 7) has a structure similar to Type 2 (resp. Type 3), but the anomalies are stronger. Types 4 and 8 show an increase of the west-to-east circulation while Type 6 has a pattern inducing a circulation in a direction opposite to the climatological west to east circulation.

Table 1 Frequency of days classified as each synoptic type for JJA and for each month from June to August

As the threshold i c decreases with increasing class rank, a greater within-type variability is found. Except for the last class, which groups together the non-classifiable days, the within-type variability (Fig. 3a) is always lower than the climatological standard deviation shown in Fig. 4. This demonstrates the added value of the 8-types CTC, compared to a study based only on climatology (i.e. where all the days are ’classified’ in one class).

Fig. 3
figure 3

The Z 500 standard deviation over the days contained in the eight classes (i.e. the within-type variability). The Z 500 standard deviation for the whole sample is shown in Fig. 4. In brackets, the mean of the within-type variability over all grid points. The solid lines (Z 500 < 5,500 m) and the dashed lines (Z 500 ≥ 5,500 m) plot the Z 500 type-average. Units are metres and the 500 hPa geopotential heights are taken from a the ECMWF reanalysis and b the NCEP/NCAR reanalysis

Fig. 4
figure 4

LeftZ500 average over the sample 4,784 days from both ECMWF and NCEP/NCAR reanalyses. Right The Z500 standard deviation for the whole sample

The next section discusses the sensitivity to different parameters and approaches used in this CTC.

3.5 Sensibility tests

3.5.1 Reanalysis and spatial resolution

Figures 2b and 3b present the CTC with the same parameter settings as the reference CTC (see Fig. 2a) but with 500 hPa geopotential heights from the NCEP–NCAR reanalysis instead of the ECMWF reanalysis. Both CTCs compare very well, even though the datasets used are different and not available at the same horizontal resolution (Fig. 4). The last synoptic circulation type (including the unclassifiable days) is different but accounts for only 1% of the days. For the other classes, the number of days per class are in agreement, although there are somewhat less (respectively more) days showing a negative [Types 2 and 5] (resp. positive [Types 3 and 7]) anomaly in the NCEP–NCAR reanalysis.

Finally, the spatial resolution used in the projection (150 or 200 km instead of 100 km) as well as the use of the 700 hPa geopotential height (as Mote 1998b) instead of the 500 hPa one do not significantly impact the results of the CTC (not shown here).

3.5.2 Number of classes

As shown in Fig. 5, the skill score S is an increasing function of the number of classes K. At K = 8, there is a break in the slope of this function and for higher values of K, S slowly tends towards a value of 0.6. In this paper, we have chosen K = 8 (S = 0.417) because the 8 types-based CTC (called hereafter CTC-8) includes all the main circulation patterns. For K higher than 8, S increases more slowly and the new synoptic circulation types found for K higher than 8 represent only a few per cent of the days.

Fig. 5
figure 5

The skill score S as a function of the number K of classes by using the ECMWF data set (in blue) and the NCEP–NCAR data set (in red). The skill score S is computed solely over the GrIS shown in blue in Fig. 2

The CTC results for K = 15 (called hereafter CTC-15) are illustrated in Fig. 6. The first nine synoptic circulation types of the CTC-15 (accounting for 92%) compare very well with the types found for K = 8 (see Fig. 2a). Only Types 10–12 are significantly new in the CTC-15 but they account for only 6% of the days. The other types (13–15) have structures similar to the 8 types found in CTC-8. There are more days in the first three classes for K = 8 than for K = 15, because part of these days are classified in other classes in CTC-15. For example, the Types 3 + 7 (resp. Types 2 + 5) of the CTC-8 are split into Types 3 + 7 + 12 + 15 (resp. Types 2 + 5 + 13) of the CTC-15. Finally, the days that ressemble to the climatology (contained in Type 1 of CTC-8) are shared in Type 1 and in Type 4 of CTC-15.

Fig. 6
figure 6

The same as Figs. 2a and 3a but for a number K of classes fixed to 15. The CTC is built using the ECMWF analysis

In agreement with Figs. 5 and 6b clearly shows that the standard deviation is lower in the classes for K = 15 (S = 0.473) than for K = 8 (see Fig. 3a).

3.5.3 Integration domain

Table 2 lists the skill scores for the CTCs using the four integration domains plotted in Fig. 7. These four CTCs differ only by the spatial dimensions of their integration domain. The 100 km-regular grid, the projection as well as the area for which the skill score S is computed (shown in red in Fig. 7) are the same. The boundary of Dom2 is one grid point away from Dom1 boundary. The Dom3 (resp. Dom4) boundaries are two (resp. three) grid points away from Dom1 boundary. Note that Dom1 is the domain used in the reference CTC presented in Sect. 3.4.

Table 2 The skill score S as a function of the number K of classes and the integration domain drawn in Fig. 7
Fig. 7
figure 7

The boundaries of four different integration domains drawn on a latitude/longitude grid. The grey scale shows the land/sea mask used in the ECMWF reanalysis. The GrIS mask used to compute the skill score S for the four domains is in red

On the one hand, as shown Fig. 8a (for Dom4), using a larger domain than Dom1 does not significantly change the results of the 8-Types based CTC. The first seven circulation types in the reference CTC are similar to these one presented in Fig. 8. Differences in the construction of classes occur because the similarity between 2 days is computed on a larger domain in Dom4 than in Dom1.

Fig. 8
figure 8

The same as Figs. 2a and 3a but using Dom4 (shown in Fig. 7) instead of Dom1. The CTC is built using the ECMWF analysis

On the other hand, the skill scores (listed in Table 2) obtained with Dom4 are lower than those obtained with Dom3, and so on down to Dom1. A lower similarity in the classes for Dom4 than Dom1 is confirmed in Fig. 8b, which shows higher standard deviations than in Fig. 8a. This is because a large integration domain takes more grid points into account in the similarity index computation outside the zone of interest.

3.5.4 Similarity index

Lund (1963) was first to design an automatic CTC based on a similarity index. His similarity index is based on the spatial correlation between two Z 500 surfaces, i.e.

$$ i_{\text{lund}}(d_1,d_2) = { \frac{\sum_{i,j=1}^{N,M} ( Z_{500}(i,j,d_1) - \overline{Z_{500}(d_1)} ) (Z_{500}(i,j,d_2) - \overline{Z_{500}(d_2)} )}{\sqrt{\sum_{i,j=1}^{N,M} ( Z_{500}(i,j,d_1) - \overline{Z_{500}(d_1)} )^2 \sum_{i,j=1}^{N,M} ( Z_{500}(i,j,d_2) - \overline{Z_{500}(d_2)} )^2 }}} $$
(8)

where \(\overline{Z_{500}(d_1)}\) is the Z 500 average of day d 1 over all grid points. The Kirchhoffer score (Kirchhofer 1973) is a derivative of the correlation coefficient introduced by Lund (1963).

Using the value of 2 for the constant c of i c in the reference CTC is justified in Table 3 noting that the exact choice of c within the listed range does not significantly affect the CTC results. The results obtained with i lund are worse in our case because the correlation coefficient compares only the Z 500 surfaces pattern and not their height (see 2). Consequently, two parallel Z 500 surfaces are classified as similar with i lund unlike our index. However, the (Lund 1963) approach gives better results if applied to all the days of a year; our approach will classify the Z 500 surface by season because the 500 hPa geopotential height is very different in summer compared to winter.

Table 3 The skill score S as a function of the number K of classes and the similarity index

Finally, the comparison presented in Table 4 justifies the use of a similarity index decreasing with the class rank rather than a fixed index i c for all classes. To classify all the days with a fixed similarity index, the algorithm has to take a lower index for the first classes compared to the one used with the decreasing index approach. Therefore, the within-class standard deviation is generally greater with the fixed index approach. In addition, almost all the days are categorized in the first three synoptic types (accounting for 95% of the sample) and the other classes (accounting for about 1%) contain only the exceptional days.

Table 4 Percentage of elements, mean within-class standard deviation and index for each class of the reference CTC presented in Fig. 2

It appears in Table 4 that the similarity index used to build the last class is zero. Indeed, to avoid a (too) low similarity index threshold when classifying all the days, our algorithm bypasses the similarity threshold for the last class if the number of unclassified days, after the k classes are built, is lower than 0.5% of the sample (i.e. about 1 day for two summers). The exact choice of this threshold does not significantly affect our CTC results.

4 Atmospheric circulation and surface melt

4.1 Surface responses to mid-tropospheric circulation changes

Figure 9 shows the anomalies of daily mean 3m-tempeature and daily meltwater production simulated by the MAR model for each class of Fig. 2. Figure 9 emphasizes the strong relationship between surface conditions simulated by MAR above the ice sheet and the circulation types of the reference CTC shown in Fig. 2. Types 3 and 7 are characterised by positive Z 500 anomalies, and are associated with warm events at the surface, while conditions colder than normal occur on the days categorized as Types 2, 5 and 8. The anomaly patterns of Types 4 and 6 are reproduced in the 3m-temperature anomalies. The 3m-temperature anomalies are often smaller along the ice sheet margin (e.g. in Type 7) than in the percolation/dry snow zone, because the surface temperature of melting snow (or refreezing snow) in the ablation zone is limited to the melting point. In most of the synoptic types, the surface is more sensitive to circulation changes along the western coast (e.g. in Type 5) than along the eastern coast. Indeed, the general west-to-east circulation produces an onshore flow (advecting cold or warm air masses following the flow direction) along the western margin of the GrIS, while there is mostly subsidence (favouring only cold air drainage from the summit of the ice sheet) along the eastern GrIS margin.

Fig. 9
figure 9

The same as Fig. 2 but for a the daily mean 3m-tempeature and b the daily meltwater production simulated by the MAR model (about 50% of this meltwater refreezes and does not reach the ocean). The daily means are averaged over the days categorized into each synoptic type of the reference CTC (see Fig. 2) and are compared with the average of all days. In brackets, the average of the anomaly (in red) and the within-type variability (in green) over the GrIS grid points. The MAR results are presented on a 25 km-regular grid as in (Fettweis 2007)

Surface melt is highly correlated with near-surface temperature, as shown in Fig. 10. The regional response of meltwater production to atmospheric circulation changes is weaker, because surface melt is only affected by surface temperature changes near the melting point isotherm (i.e. in the ablation zone along the ice sheet margin).

Fig. 10
figure 10

Time series of the mean GrIS 500 hPa geopotential height and 2m-temperature from both ECMWF and NCEP–NCAR reanalyses as well as the daily 3m-temperature, solar incoming radiation, meltwater production and meltwater run-off simulated by MAR averaged over the 1958–2009 period. The time series are centred (i.e. with a mean of 0) and normalised (i.e. with a standard deviation of 1) to facilitate comparison

The within-class standard deviation for 3m-temperature and meltwater is of the order of the standard deviation of the whole sample (in contrast to Z 500 for Types 1–7 as shown in Fig. 2), because the surface conditions can be very different for similar atmospheric circulations. On the one hand, the thermal inertia of the snow (e.g. the snow temperature has to reach 273.15 K before melting starts) offsets the surface response to atmospheric changes. On the other hand, as listed in Table 5, a mid-tropospheric circulation anomaly at the start of the melt season, when solar radiation is maximum but albedo is high because of the long-lasting winter snowpack, has a different impact than the same anomaly at the end of the ablation season, when low albedo zones (bare ice) appear and the incoming solar flux is decreasing.

Table 5 The total GrIS available meltwater (in km3) in June, July, August and the total in JJA for the days categorized in the eight synoptic types as well as for all the summer days of the 1958–2009 period

Phase shifts between the different times series are illustrated in Fig. 10:

  • Near-surface temperatures (both from the MAR model and ECMWF reanalysis) are in phase with meltwater production although, at the end of August, surface melt still occurs in the lower regions of the ablation zone, while winter conditions return to other parts of the ice sheet.

  • There is a lag of about 15 days between meltwater production and run-off. Meltwater run-off occurs in the MAR model when an internal and accumulated surface meltwater threshold is reached according to (Lefebre et al. 2003). In average, only about 50% of the available meltwater reaches the ocean in MAR. Despite this phase shift between the time series, the interannual run-off variability is almost fully explained (r = 0.98) by the cumulated meltwater production.

  • Both reanalyses show the same daily variability of Z 500 superimposed on the sinusoidal seasonal cycle. The near-surface temperature and meltwater time series are somewhat affected by the Z 500 variability, but the run-off time series is fully smoothed. Finally, there is a clear phase shift between Z 500 and the near-surface temperature. This explains the higher within-class standard deviation for 3m-temperature and meltwater compared to Z 500.

Hereafter, we will call a synoptic type warm (resp. cold) if this type induces on average warmer (resp. colder) surface conditions than normal, according to Fig. 9.

4.2 Interannual variability

The time series of the main SMB components as well as the near-surface temperature simulated by MAR are plotted in Fig. 11. In agreement with Hanna et al. (2008) and Rignot et al. (2008), the GrIS gains mass at its surface from the beginning of the 1970s until the middle of the 1990s, due to a conjunction of cold and wet years, and loses mass in the dry and warm 1960s and in the very warm 2000s summers. Surface melt has been higher than normal in the early 1960s and over the last 10 years; the current abnormal melt is very likely due to anthropogenic warming i.e. the human activities as concluded by Fettweis (2007) and Hanna et al. (2008). The PDD time series clearly show the highest correlation with the surface melt time series, while the correlation of surface melt with the mean JJA surface temperature is already very good.

Fig. 11
figure 11

Top Time series of main annual ice sheet mass budget components (SMB ∼ snowfall - run-off) simulated by MAR. Units are km3yr−1. Bottom Time series of the GrIS JJA mean 3m-temperature, cumulated Positive Degree Days (PDD) and meltwater production. The mean JJA Z500 simulated by ECMWF reanalysis is also shown. Correlations with run-off (mauve) and meltwater production (black) are listed in brackets. The time series are centred and normalised, and the 5-year running mean for meltwater is plotted as a dashed black line

Since 1998, the 500 hPa-general circulation seems to favour warm events during the summers, as shown in Fig. 12 and in Table 6. Both reanalyses show a decreasing number of cold synoptic types. Before the 2000s, interannual variability is too high to detect trends. The summers with the highest number of days classified in the warm types are 2007 and 2008 (followed by 1960 and 1998). This explains why a record surface melt is simulated by MAR for the two last summers.

Fig. 12
figure 12

Top Time series of the number of days where the mid-tropospheric circulation promotes cold events (Type 5 from the CTC of the Fig. 2a in dark blue), moderate cold events (Type 2 + 8, in light blue), moderate warm events (Type 3, in light red) and warm events (Type 7, in dark red). The days classified as Type 1 (similar to the 1958–2009 climatology) are not taken into account. The time series are centred. Dashed lines represent the 5-year running mean of the unified red (resp. blue) times series. Below The same but with the CTC of Fig. 2b built using the NCEP–NCAR reanalysis

Table 6 Number of days in each class of the reference CTC (see Fig. 2a) for each summer (JJA) of the 1958–2009 period

Minima in the MAR surface melt time series occur in 1983 and 1992 after the respective eruptions of El Chichon and Pinatubo, as found by Hanna et al. (2005) and Fettweis (2007). By injecting large amounts of aerosols in the atmosphere, volcanic eruptions reduce the amount of solar energy reaching the Earth surface, lowering the surface temperature. On the other hand, the summer with the highest number of days classified in the very cold types is 1983. The record low JJA 3m-temperature found in 1983 can thus be explained by a combination of the El Chichon eruption and of an abnormal high number of circulation-induced cold events. In 1992, the impact of the mid-tropospheric circulation is less significant.

The North Atlantic Oscillation (NAO) represents the dominant mode of regional atmospheric variability around Greenland (see Rogers 1997; Bromwich et al. 1999) and is quantified by the NAO index. Typically, the NAO index is based on the monthly normalised pressure difference between Gibraltar and Reykjavik. However, we here use the NAO index from the Climate Prediction Center (CPC) because it is based, as our CTC, on the daily (00Z) 500 hPa height anomalies (see http://www.cpc.ncep.noaa.gov/ for more details). During a negative NAO phase, the weaker Icelandic low favours (southerly) warm air advection along the southwest coast of Greenland and over the GrIS. Although the GrIS sensitivity to the oscillation is strongest in winter (Fettweis 2007), we nevertheless found positive (resp. negative) correlations between the mean JJA NAO index and the number of summer days with an atmospheric circulation that induces cold (resp. warm) events in agreement with Hanna and Cappelen (2003). This correlation is 0.78 (resp. −0.8) with the number of summer days of Type 2 or Type 5 (resp. Type 3 or Type 7). While it is questionable to directly link the NAO index to the number of days in synoptic classes, these correlations show that the interannual variability found in our CTC is linked to the NAO.

4.3 Daily variability

Figure 13 shows that, in both reanalyses, the daily melt and near-surface temperature anomalies within a single season are highly correlated with 500 hPa circulation variability although the surface responses to synoptic variability can be very different due to the surface inertia. In addition, it should be noted that

  • the warm or cold nature attributed to a synoptic type is computed by comparison to the JJA average, which does not prevent a 500 hPa circulation classified as a warm (cold) synoptic type from inducing several days colder (warmer) than normal conditions at the surface.

  • there is sometimes a delay of 1 day between temperature change and circulation change because the temperatures shown in Fig. 13 are daily means and the circulation type is based on a snapshot at 12 h UT.

  • the CTC’s using the ECMWF reanalysis (see Fig. 2a) and the NCEP–NCAR reanalysis (see Fig. 2b) are slightly different as shown by the number of days classified in each synoptic type. This explains why a day can not be classified as the same type using ECMWF-based CTC or NCEP-based CTC.

Fig. 13
figure 13

Top Time series of daily 3m-temperature (in red) and meltwater (in blue) anomalies simulated by MAR during summer 2008. The anomalies refer to the 1958–2009 JJA average. The time series of the daily 1958–2009 mean 3m-temperature and meltwater (shown in Fig. 10) are also plotted. Below Time series of the synoptic types of each day from the reference CTC shown in Fig 2. Note that the synoptic types (while they are categories and not values) are rearranged over the Y-axis according to the caption order (from colder to warmer types). Plots like this Figure are available for each summer of the 1958–2009 period in the electronic supplementary material

In 2008, 3m-temperatures are higher than the JJA mean from the beginning of June until June 20. The anomaly peaks June 14, after five consecutive days with the warm Types 3 and 7. Around June 20, the synoptic type is near the climatology as well as the surface temperature. Afterwards, there is a second warm period, characterized by Types 3, 6 and 7, until the first decade of July. At the end of July, there is a third warm period lasting 15 days, with a succession of warm synoptic types, advecting warm air as far as northern Greenland. Around August 10, Type 2 brings back this mass of warm air to the south, inducing temperatures still higher than the JJA normal. Mid-August, Type 8 atmospheric circulation decreases the temperature below the climatology. Finally, at the end of August, a fourth warm period occurs over northern Greenland induced by Type 6 circulation. The 2008 summer has the highest number (15) of days classified as Type 6 (representing a warm anomaly in the North Greenland). This result is in agreement with Tedesco et al. (2008b), who detected extreme snowmelt in Northern Greenland through passive microwave data. Finally, Fig. 13 highlights the very high correlation between near-surface variability and surface melt. The surface melt response to circulation changes is nevertheless more smoothed than the 3m-temperature responses, due to snowpack inertia.

Both reanalyses (see Fig. 2) classify the four warm periods of the 2008 summer as Types 3 or 7. However, there are more Type 7 days in the NCEP-based CTC; this is because the associated warm anomaly is smaller in the NCEP-based CTC than in the ECMWF-based CTC. Consequently, Class 7 is larger in the NCEP-based CTC. While the NCEP–NCAR reanalysis is not used to force the MAR model, the MAR simulated 3m-temperature anomalies in 2007 are in agreement with the synoptic types derived from the CTC presented in Fig. 2b. There are five warm events in 2007 with a temperature about 5 K higher than the climatology (see Fig. 14). During theses events, the 500 hPa circulation is classified as Types 3 or 7, the warmest synoptic types in the CTC using the NCEP–NCAR reanalysis. Summer 2007 is characterized in both reanalyses by a very high number of Types 3 and 7 favouring advection of warm air mass over southwestern Greenland. This concurs with the surface anomalies found by Tedesco et al. (2008a) during summer 2007.

Fig. 14
figure 14

The same as Fig. 13 but for summer 2007

During summer 2003, an exceptional warm event occurred at the end of August (see Fig. 15). During the period 26–28 August 2003, the near-surface temperatures were 6K above the JJA average. When compared to the climatology for the end of August, the anomalies were 11 K above normal while surface meltwater production was seven times the average value for the end of August. This episode is classified as Type 7, the warmest synoptic type in the reference CTC. The days where the similarity index with 27th of August 2003 is higher than 0.80 are the 6th of June 1971, the 27–29th of August 2003, the 11th of August 2004 and the 15th of August 2006 (see Figs. 15 and 16). Every time such a circulation occurred, the temperature was on average about 8 K above the climatology. In addition, the Type 7 circulation pattern lasted for 6 days in 2003, something that otherwise happened only in 1966, 1971, 2007 and 2008. This confirms the exceptional nature of late August 2003.

Fig. 15
figure 15

The same as Fig. 13 but for August 2003, 2004 and 2006

Fig. 16
figure 16

The same as Fig. 13 but for summers of 1971 and of 2009

The highest melt event occurred on the 27th of June 2007 (classified as Type 3 in the reference CTC) and resulted from the persistence of circulation-induced warm events since the beginning of June. The second highest melt event (i.e. three times the JJA mean) occured on the 16th of July 2009. The 500 hPa atmospheric circulation of mid-July 2009 (see Fig. 16) are not unique but this is the only time in the period 1958–2009 that a Type 7 circulation persists during 9 days in the middle of July, when surface melt is near its seasonal peak.

4.4 Conclusion

The automatic CTC over Greenland presented here is based on the 500 hPa geopotential height, and aims at quantifying the impact of mid-tropospheric circulation variability on GrIS surface melt simulated by the MAR model. In contrast to the (Lund 1963) classification, this CTC is based on the height of the Z 500 surfaces rather than on their pattern. In addition, it includes a method for reducing within-type variability during the construction of synoptic types.

The CTC shows that the dominant mode of regional atmospheric variability around the GrIS is linked to NAO and that surface melt anomalies are highly correlated to the general circulation. It explains why a record surface melt was observed during the summers of 2007 and of 2008. The CTC also highlights the changes that have been affecting the JJA atmospheric circulation since the end of the 1990s. Indeed, the 500 hPa-general circulation seems to favour warm summer events over the ice sheet more now than before. As found by Philipp et al. (2007) for Europe, these circulation changes significantly contribute to the GrIS warming observed since the 1990s (Steffen and Box 2001; Chylek et al. 2006; Hanna et al. 2008).

An interesting future development of the present research would be to link the precipitation variability over the GrIS to the atmospheric circulation types by applying our CTC over the whole year (i.e. season by season). In addition, we plan to apply our CTC to the outputs of the atmosphere-ocean general circulation models (AOGCMs) of the Fourth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR4), to detect potential changes in the future atmospheric circulation over Greenland. The CTC should be applied to the outputs from the “Climate of the Twentieth Century Experiment” (20C3M) and compared with the CTC obtained with both ECMWF and NCEP–NCAR reanalyses to assess the AOGCM’s efficiency to reproduce the current general circulation. Finally, it should be interesting to understand why the atmospheric circulation is changing over Greenland due to outgoing global warming.