The many studies investigating the future change of the Greenland Ice Sheet surface mass balance from climate model output exhibit a wide range of projections. This study makes projections from the Coupled Model Intercomparison Project phase 3 models used in the Intergovernmental Panel on Climate Change Fourth Assessment Report and explores the underlying physical processes behind their spread. The projections are made for three Special Report on Emissions Scenarios, B1, A1B, and A2, with a focused analysis on the A1B scenario. The estimate in the study suggests that about 60% of the intermodel difference in the twenty-first-century ablation rate change under the A1B scenario is accounted for by the global annual mean temperature change. In the current study, other processes responsible for the spread in model projections are investigated after excluding this global effect. A negative correlation (−0.60) was found between the simulated summer temperature bias over the Greenland Ice Sheet under present-day conditions and the ablation rate increase during the twenty-first century, partly because maximum warming over ice is approximately limited to the melting temperature. Models with relatively larger ablation rate increases during the twenty-first century exhibit greater warming with a greater reduction in sea ice cover. The authors found that these models also simulate relatively cooler summer conditions in high latitudes with more sea ice cover in the late twentieth century, suggesting the importance of sea ice feedbacks. Also, an anticorrelation (−0.75) is found between weakening of the Atlantic meridional overturning circulation and the ablation rate increase during the twenty-first century. The relation in the model spread between the twenty-first-century ablation change and the late twentieth-century climate conditions is then used to investigate the impact of model bias on the multimodel ensemble of projections. The result suggests that the models’ underestimation of present-day sea ice concentration near the coast of Greenland may cause an underestimation of future Greenland Ice Sheet ablation rate increase in the ensemble projection. These results emphasize the importance of correctly simulating present-day conditions and understanding the underlying multiple physical processes behind the intermodel difference to reduce the uncertainty in future projections.
Future sea level rise is of great concern for society as it may cause many hazards such as coastal inundation, erosion, submergence of small islands, destruction of freshwater resources, ecosystem losses, and increasing vulnerability to other natural hazards (Nicholls et al. 2007; Mimura et al. 2007; Nicholls and Cazenave 2010). Bindoff et al. (2007) reports that the global mean sea level rise over the twentieth century was about 17 cm. For the decade of 1993–2003, about half of the observed rise was attributable to the steric component, about 25% was caused by melting of glaciers and small ice caps, and only about 14% was contributed by the Greenland and Antarctic ice sheets (Bindoff et al. 2007). Enhanced contribution from the polar ice sheets is, however, suggested for 2003–08 compared to the previous decade (Cazenave et al. 2009). The projected sea level rise for the twenty-first century is 18–59 cm with an additional 10–20 cm uncertainty due to rapid dynamical changes of the ice sheets, which amounts to 18–79 cm (Meehl et al. 2007a). There is still large uncertainty in these numbers, and individual contributions need to be continuously reexamined based on the understanding of the physical processes.
Recent observations show rapid mass loss of the Greenland Ice Sheet (GrIS) (Chen et al. 2006; Cazenave 2006; Shepherd and Wingham 2007; Wouters et al. 2008; Baur et al. 2009; Velicogna 2009; Khan et al. 2010). Huybrechts et al. (2004) list four sources of uncertainty in model estimates of future ice sheet change: 1) long-term background evolution, 2) surface mass balance (SMB) changes, 3) ice dynamic response to these mass balance changes, and 4) rapid dynamical changes in ice flow. Quantification of the first source requires the long-term integration of an ice sheet model over glacial–interglacial cycles with proper representation of isostatic adjustment (e.g., Abe-Ouchi et al. 2007). Huybrechts and de Wolde (1999), van de Wal et al. (2001), and Huybrechts et al. (2004) concluded that the contribution of the third source is not negligible for the GrIS, but is secondary by the end of the twenty-first century. The last source refers to ice flow acceleration caused by processes poorly represented in the current ice sheet models, such as inland transmission of stress after ice front thinning or breakup (Alley et al. 2005; Howat et al. 2007; Pritchard et al. 2009; Nick et al. 2009) and enhanced basal lubrication after surface melting (Zwally et al. 2002; van de Wal et al. 2008; Joughin et al. 2008; Shepherd et al. 2009; Bartholomew et al. 2010). In the current study, we focus on the uncertainty in the second source, that is, the future change of the GrIS surface mass balance.
Observations suggest that there has been an increase in GrIS melting extent over time (Mernild et al. 2009), and studies using regional climate models suggest a decrease of GrIS SMB in recent decades (Fettweis 2007; Hanna et al. 2008; Rignot et al. 2008; van den Broeke et al. 2009). Previous studies have shown a strong dependency of the future SMB projection on the given climate model. Ohmura et al. (1996) and Wild and Ohumra (2000) performed time-slice experiments using two T106 (at that time, relatively high resolution) atmospheric general circulation models (AGCMs) with boundary conditions taken from a CO2-doubling transient experiment using a low-resolution (T21) atmosphere–ocean GCM (AOGCM). These two different models produced distinctly different SMB estimates. As the narrow melting zone is not sufficiently resolved, even with T106 AGCMs, Wild et al. (2003) interpolated the temperature anomaly under doubled-CO2 concentration onto a 2-km horizontal grid and added it to the observed climatology so as to estimate ablation more accurately. Suzuki et al. (2005) applied a similar method for two transient scenario simulations. The contribution of the GrIS to sea level change is positive in Suzuki et al. and negative in Wild et al.. To compute the SMB for a large number of scenarios, Bugnion and Stone (2002) used a snowpack mass balance model on a 20-km horizontal grid driven by output from an earth system model of intermediate complexity (EMIC). In their study, the sign of the GrIS contribution to sea level change depends on the idealized scenario. Church et al. (2001), Huybrechts et al. (2004), Gregory and Huybrechts (2006), and Meehl et al. (2007a) used a pattern scaling technique. In this technique, spatial patterns of temperature and precipitation changes simulated by high-resolution A(-O)GCMs are scaled by the annual mean ice sheet average or global average temperature and precipitation changes simulated by multiple lower resolution AOGCMs or EMICs that are run continuously under various scenarios. In these studies, the SMB was evaluated on a 20-km horizontal grid and the best available information on climate change is incorporated, although the number of the high-resolution models is limited. Fettweis et al. (2008) applied a regression technique to the multi-GCM output. The regression equation is derived from a regional climate model and predicts the SMB as a function of temperature and precipitation changes in specific regions of the GrIS.
In most previous studies, the SMB was estimated diagnostically from climate model output. While it exhibits a wide range of projections and the scenario dependency was investigated to some extent, the underlying physical processes behind the model spread have not been identified. In the current study, we attempt to reveal the processes relevant to the range of GrIS SMB projections. Although assessing the ice sheet contribution to sea level change more precisely requires the inclusion of ice dynamics, this is not targeted in this study. However, an important aim of this study is to provide information on how the results of such extended studies using the dynamic ice sheet models (Charbit et al. 2008; Vizcáino et al. 2010; Ridley et al. 2010; Graversen et al. 2011) could be affected by various processes through uncertainties in the SMB change.
Once important processes are identified, we aim to make use of this information to interpret the ensemble of future projections. In general, a multimodel ensemble of future projections may be evaluated based on the performance of the present-day simulation (e.g., Franco et al. 2010). In some studies, a weighting was assigned for each model to obtain the best estimate of the projections and to reduce the range of uncertainty (e.g., Murphy et al. 2004; Schmittner et al. 2005). In such an approach, an intuitive assumption is made that models that reproduce the current conditions well are expected to project the future well. Alternately in this study, we first identify important variables and their spatial patterns in the present-day simulation that are related to the intermodel difference in future projections. These spatial patterns are then linked to the model bias in the present-day simulation to assist the interpretation of the multimodel ensemble of future projections.
This article is organized as follows. In the next section, the multimodel dataset, observational data, and the method used in calculating the SMB are described. In section 3, the projected SMB changes are characterized, and physical processes relevant to the spread of the projections are identified. Future ensemble projections are then evaluated in terms of model bias in the present-day simulation. The remaining uncertainties are discussed in section 4, followed by conclusions in section 5.
2. Data and method
The GrIS surface mass balance (SMB) is estimated from the output of 25 GCMs. The twentieth-century simulation and future projections based on three Special Report on Emissions Scenarios (SRES) (Nakićenović et al. 2000) are taken from the World Climate Research Programme (WCRP) Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel dataset (Meehl et al. 2007b). The three scenarios—B1, A1B, and A2—correspond to CO2 concentrations of about 550, 700, and 820 ppm in the year 2100, and some models were integrated continuously for another 100–200 years with fixed atmospheric constituents.
The SMB is only calculated over the current topography of the GrIS without allowing for changes in elevation through time. The ablation rate is calculated as a function of summer (June–August) surface air temperature using the so-called temperature index (TI) method according to Ohmura et al. (1996) and Wild et al. (2003) in which ablation increases linearly with temperature above −1.8°C. It contains contributions from both surface melting and refreezing. A physical basis for this formulation is given by Ohmura (2001). While the TI method is also used in Kiilsholm et al. (2003) and Suzuki et al. (2005), many other studies have used the positive degree-day (PDD) method (e.g., Gregory and Huybrechts 2006; Meehl et al. 2007a), and a discussion of these schemes is given in section 4. Summer temperature is obtained by adding observation-based gridded temperatures in the fine resolution 0.02° latitude × 0.05° longitude National Survey and Cadastre (KMS) grid (Ekholm 1996) in a reference period and model-simulated temperature changes from that reference period. This kind of anomaly method has been widely used in many studies on ice sheet mass balance (e.g., Gregory and Huybrechts 2006; Meehl et al. 2007a; Vizcáino et al. 2008; Ridley et al. 2010). The observational data are taken from Calanca et al. (2000), and model data are interpolated to the KMS grid. This resolution (~2 km) is much higher than that used in typical previous studies (~20 km). While Calanca et al. notes the representative period of 1951–60, we take the longer period, 1901–60, as the reference period for the model data so as to reduce the stochastic effect of internal variability. The results are confirmed as insensitive to this particular choice of reference period because the variability in the twentieth century is much smaller than the change in the twenty-first century (not shown). The present-day melting rate computed internally in coarse-resolution GCMs tends to be overestimated because the ice sheet is smoother and has a lower elevation compared to that in reality (Thompson and Pollard 1997; Glover 1999; Kiilsholm et al. 2003). On the other hand, the anomaly method taken here could, in principle, allow inclusion of the first-order effect of the topography of the ice sheet. This also minimizes the direct influence of GCM biases in estimating ablation. While Gregory and Huybrechts (2006) reduced the effect of coarse resolution of the CMIP3 GCMs by using spatial patterns simulated with only high-resolution GCMs, we reduce the effect by using observed background climatology on very high resolution.
The accumulation rate is calculated by simply interpolating the GCM-simulated snowfall rate on the same KMS grid. Thus, sublimation is ignored in the current study, which may result in an overestimation of the accumulation by about 10% (Wild et al. 2003). The smoothed topography in GCMs likely overestimates precipitation inland of the ice sheet since there is an insufficient barrier against moisture transport at the margins (Wild et al. 2003), but a fraction of solid precipitation may be underestimated as more precipitation is diagnosed as rain rather than snow where the surface elevation is underrepresented by smoothing. Improvement of the accumulation estimate is in progress (Bales et al. 2009; Burgess et al. 2010), and Ettema et al. (2009) also emphasizes the importance of resolution in correctly simulating the accumulation. These factors introduce uncertainties into the absolute magnitude of the projected SMB change. Nevertheless, our essential conclusion does not change because ablation change is dominant, as suggested in previous studies (e.g., Ridley et al. 2005; Vizcáino et al. 2010), and thus is the focus of the current study. Note that virtually the same procedures for calculation of ablation and accumulation were followed by Suzuki et al. (2005) in estimating the SMB of the Greenland and Antarctic ice sheets for two of the CMIP3 models.
a. Surface mass balance projections
Figure 1a shows projected GrIS SMB changes represented in the equivalent rate of sea level rise for the A1B scenario with respect to the 1900–49 average. The intermodel differences become pronounced by 2100, and the range of the spread remains at a similar level thereafter due to the stabilization scenario. By 2300 the rate differs by as much as a factor of 4 between the maximum (ECHAM5/MPI-OM) and minimum (PCM) model results. (see footnote in Table 4) The contributions of snowfall and ablation to the SMB change are shown in Figs. 1b and 1c, respectively. It is clear that ablation dominates over snowfall in determining both the sign and spread of the net change. The 5%–95% range of changes in the equivalent rate of sea level rise between 1980–99 and 2090–99 is 0.2–2.0 mm yr−1, comparable to 0.3–1.9 mm yr−1 estimated in Meehl et al. (2007a).
Information on the SMB change alone does not allow estimation of the contribution to sea level change. Following the approach of Suzuki et al. (2005), we assumed that the total mass of the ice sheet is in balance for the 1900–49 average (cf. Huybrechts et al. 2004) and that there is no contribution to sea level change thereafter from non-SMB components. Contribution to sea level change is then obtained by integrating the SMB rate change from 1900 (Fig. 2). The results are relatively insensitive to the choice of averaging period (not shown) as there was little SMB change in the twentieth century. The 5%–95% range of the GrIS sea level contribution between 1980–99 and 2090–99 for the A1B scenario is 1.9–12.9 cm, which is larger than the 1–8 cm stated in Meehl et al. (2007a). The noticeable difference results from a small difference in the total SMB rate change and the speed of the rate changes over the integration period. Intermodel spread in the sea level contribution grows substantially after 2100 as it is the integrated rate of change. As for the 13 models available until 2300 for the A1B scenario, the 5%–95% range of the GrIS sea level contribution with respect to the 1980–99 average is 3.4–12.7 cm at the end of the twenty-first century and becomes 17.5–70.5 cm by the end of the twenty-third century. A comparison with Meehl et al. for the three scenarios is summarized in Table 1, and time and scenario dependence are summarized in Table 2. Hereafter, our investigation focuses on the A1B scenario, in which the available number of models and the projected range of the uncertainty is large (Table 2). In addition, we restrict our analysis to the change in the twenty-first century (from 1980–99 to 2080–99), when the steepest increase of the ablation rate occurs.
b. Relation to global and Greenland temperature changes
As a first step in exploring reasons for the large model difference, correlations r between variables of interest are computed. The correlation matrix for the twenty-first-century changes is presented in Table 3. Consistent with the dominance of ablation over accumulation in Fig. 1, there is a very high correlation between the equivalent rate of sea level rise (Δnet) and that due to ablation alone (Δabl). As expected from the empirical formula that we applied, there is a high correlation between the ablation rate (Δabl) and GrIS summer temperature . The nonperfect correlation is explained by the spatial heterogeneity of temperature change over the GrIS. The correlation between GrIS and global annual mean temperature change is 0.67, supporting the idea that some local processes, such as sea ice albedo feedback, also contribute to the temperature change near Greenland. While there is an almost perfect correlation between summer and annual-mean global average temperature changes and , the correlation between summer and annual GrIS temperature changes and is about 0.74. This also implies that there are some local processes in determining the seasonal change near Greenland. As a result, the global annual-mean temperature change accounts for about 60% (r = 0.78) of the ablation rate change (Fig. 3a). The global annual-mean temperature change is controlled by climate sensitivity and the rate of ocean heat uptake, and these variables have been investigated extensively in previous studies (e.g., Soden et al. 2008; Yokohata et al. 2008). Thus, we focus on other factors that affect the change in the ablation rate.
c. Relation to GrIS temperature bias
We now turn our attention to the local GrIS temperature change. Figure 4a shows a 24-model ensemble mean of projected summer temperature changes during the twenty-first century for the A1B scenario. It features amplified warming in central Greenland, and anomalous isothermal contours loosely follow the terrain contours.
In typical increased CO2 experiments, warming is confined to lower altitudes in a strongly stratified high-latitude atmosphere (Manabe and Wetherald 1975). Also, amplified warming caused by sea ice albedo feedback is expected to have a larger impact on the neighboring coast rather than far inland. Therefore, neither stratification nor albedo change account for enhanced warming in the high-elevation area seen in the ensemble mean. It is likely that the suppression of warming in the low-elevation area is due to increased land snow/ice melting and sea ice melting in GCMs, as pointed out by Gregory and Huybrechts (2006). The increased available energy during the twenty-first century is consumed as latent heat by snow/ice melting, and the surface temperature remains near the melting point in those regions. The ablation rate increase diagnosed directly from the GCM-simulated temperature without using the anomaly method supports this interpretation (Fig. 4b). This area of increased ablation is much broader than that estimated using the anomaly method (Fig. 4c). Note that larger and smaller ablation rate increases in the interior and coastal regions, respectively, in Fig. 4b compared to Fig. 4c primarily reflect the initial temperature bias in GCMs (Fig. 4d). These figures underscore the fact that the warming simulated by GCMs includes the result of snow/ice melting computed internally in GCMs and raises concerns about the use of GCM temperature as input to ablation calculation. The use of empirical formulas in regions of melting in both GCMs and the fine KMS grid is, however, physically consistent because the formula fits the observed relation between temperature and ablation over the melting surface. Its use in regions of melting in GCMs, but initially nonmelting in the KMS grid because of higher elevation, may result in an underestimation of potential ablation. Consequently, the SMB projection using the anomaly method is still subject to model bias in the initial GrIS temperature. Therefore, this problem has occurred in many previous studies, and it is one of the key points that we highlight below.
Figures 5a and 5b show summer surface air temperature changes in two groups of model results having larger (>60 mm yr−1, L group) and smaller (<−60 mm yr−1, S group) ablation rate increases, respectively, with respect to the regression line in Fig. 3a. The threshold was determined so that each composite retains five members. Note that three models (FGOALS-g1.0, PCM, and GISS-EH) are not considered as candidates for the two groups because either sea ice data are unavailable or sea ice distribution in the late twentieth century is too unrealistic for later analysis. The L group exhibits greater amplification of warming over Greenland than the S group. In addition, more enhanced warming and ablation rate increases are seen in northern Greenland in the L group compared to the S group (not shown). Figures 5c and 5d show composites of GCM-simulated summer temperature (before applying the anomaly method) in the late twentieth century for the L group and the S group, respectively. They show that the models with large increases in the ablation rate actually have a colder climate than the models with small increases in the ablation rate. This is because models with a colder climate in the late twentieth century require more warming to occur until the surface temperature reaches the melting point. To quantify the effect of this model bias, a correlation is computed between the GCM-simulated GrIS summer surface air temperature bias from the data of Calanca et al. (2000) and the ablation rate increase normalized by the global annual mean temperature change (normalized ablation rate change) during the twenty-first century (Table 4; Fig. 3b). The correlation coefficient is −0.60 and is statistically significant at the 1% level, suggesting that GCM biases in Greenland temperature affect the degree of warming and projected ablation change. The fact that the increase in ablation rate diagnosed directly from the GCM-simulated temperature, without using the anomaly method, is larger in the S group than in the L group (not shown) supports the interpretation that more latent heat is consumed in melting and potential warming is suppressed in the S group GCMs.
d. Relation to ocean circulation change
Here we discuss the influence of warming outside Greenland on the temperature change over the ice sheet. A composite difference in summer temperature change during the twenty-first century between the L group and the S group is shown in Fig. 6. Warming is not confined to the ice sheet, but rather shows a large-scale feature in the northern North Atlantic Ocean. This suggests that the temperature change around Greenland and to the south of Greenland may be important in explaining the model difference in ablation change.
The reason for strong differences in temperature change to the south of Greenland between the composites is more easily understood when composites for the L group and the S group are presented separately (Figs. 7a and 7b). The S group actually shows cooling to the south of Greenland. This region of cooling under increased CO2 is a feature associated with the weakening of the Atlantic meridional overturning circulation (AMOC), which presently transports heat from low to high latitudes. To verify this, we made composites for models that exhibit weakening of the AMOC by less than 3.5 Sv (1 Sv ≡ 106 m3 s−1) and by more than 5 Sv, as shown in Figs. 7c and 7d (see also Table 5). These thresholds were determined so that each composite retains initially five members. We exclude two models (INM-CM3.0 and FGOALS-g1.0) because they exhibit an unrealistic AMOC structure in the late twentieth century and their changes in the streamfunction are not interpreted as “AMOC weakening.” A very similar region of cooling is seen in both temperature change in the S group (Fig. 7b) and the group with large AMOC weakening (Fig. 7d). In addition, the L group and S group models are part of the selected models of small and large AMOC weakening, respectively. This suggests a high possibility that AMOC changes influence the GrIS surface mass balance. The largest negative temperature anomaly to the south of Greenland after subtracting the zonal mean temperature change is −2.2° and −3.1°C in the L group (Fig. 7a) and the S group (Fig. 7b), respectively, and −1.6° and −3.4°C in the group of small and large AMOC weakenings (Figs. 7c and 7d), respectively. There is a strong anticorrelation (−0.75) between AMOC weakening and the normalized ablation rate change (Fig. 3c). The normalized ablation rate increase is not, however, particularly enhanced over southern Greenland where the largest influence of AMOC weakening is expected (not shown). This is probably because the ablation rate increase is also affected by the model bias discussed in section 3c. Note that Gregory et al. (2005) reported that models exhibit AMOC weakening under a quadrupling of CO2 (4 × CO2), but not cooling. They also found a positive correlation between the control AMOC strength and the degree of weakening under 4 × CO2, but we did not find such a relation between the late twentieth-century AMOC strength and the degree of weakening during the twenty-first century under the A1B scenario (Table 5). The existence of cross-model correlations does not immediately imply causality, but the physical processes are well understood. For example, the impact of AMOC weakening on the ice sheet mass balance is well demonstrated in Vizcáino et al. (2008) in which the tripling-of-CO2 experiment (3 × CO2) exhibits much larger ablation than the 4 × CO2 experiment. Also, the effect of meltwater increase on AMOC weakening would produce a correlation of opposite sign, if any.
e. Relation to sea ice distribution
The large area of warming along the coast of Greenland in the composite difference between the L group and S group suggests that sea ice distribution may also play an important role in ablation rate changes (Fig. 6a). Indeed, a larger reduction in summer sea ice cover near the coast of Greenland during the twenty-first century is seen in the L group compared to the S-group (Fig. 8a). In contrast, a difference in sea ice thickness change between the L group and S group mainly occurs in the Arctic Ocean, rather than near Greenland (Fig. 8b). The different signs in the ice concentration difference (Fig. 8a) and ice thickness difference (Fig. 8b) within the Arctic Ocean reflect the fact that the concentration changes little with small melting when the initial ice amount is large, while the thickness changes proportionally to the melting. There is a positive correlation between the reduction in summer sea ice in the area 45°–85°N, 90°–30°W and the normalized ablation rate change (Fig. 3d). The correlation coefficient is, however, sensitive to the specification of the region. As causality is not demonstrated from this analysis, the possibility remains that other processes responsible for polar amplification of warming, such as meridional heat transport change or lapse rate feedback (Yoshimori et al. 2009), cause the warming both outside and within Greenland and create the apparent correlation.
Figure 6b shows the composite difference in GCM-simulated summer temperature (before applying the anomaly method) between the L group and S group in the late twentieth century. It clearly shows that models with a colder climate in a large area near Greenland in the late twentieth century experience greater warming during the twenty-first century. It is important to stress that the colder conditions of the late twentieth century are not confined to the ice sheet. Sea ice concentration and thickness in the L group consistently indicate colder conditions, that is, larger sea ice extent and thicker ice in the late twentieth century (Figs. 8c and 8d). Sea ice plays an important role in determining the energy budget through reflecting solar radiation with high albedo and limiting air–sea heat exchange. As these feedbacks only operate in regions where ice initially exists, it is physically plausible that models with greater ice extent near Greenland in the late twentieth century experience greater warming when the ice is removed. A qualitatively similar suggestion was made by Holland and Bitz (2003) and Boé et al. (2010). Nevertheless, it is possible that the link to the cold bias outside Greenland is apparent due to the concurrent cold bias over the GrIS, as discussed in section 3c. Qualitative aspects of the results do not, however, change even if models with large biases (>3°C) in summer GrIS temperature are excluded. Note also that there is no significant correlation (at the 10% level) between AMOC weakening and the reduction of summer sea ice cover. Nevertheless, two extreme examples (models 2 and 24) show both stronger (weaker) AMOC reduction and weaker (stronger) sea ice reduction, and thus the two variables may not be totally independent.
f. Effect of model bias on the ensemble projection
In the previous subsection, it was shown that the model-predicted future change in ablation rate over Greenland is related to the simulated late twentieth-century climate conditions over a large area of the Northern Hemisphere. Such a relation enables us to evaluate the ensemble projection in light of the performance of models in simulating the observed climate conditions. To make use of the spatial information in the model spread, we apply singular value decomposition (SVD) analysis in spatio-model space for two variables: ablation rate and summer surface air temperature changes. The SVD extracts leading spatial patterns that share the variance for the two variables across the models. As these statistical methods have been frequently used in literature, we refer readers to Bretherton et al. (1992) and von Storch and Zwiers (1999, 293–333) for their mathematical description. In addition, a mathematically identical procedure to that used here is described in great detail in Shiogama et al. (2011). The two variables are first normalized by the global annual-mean surface air temperature change to exclude the effect of the climate sensitivity difference, and then further normalized by the standard deviation of each variable averaged over the domain of the analysis to equalize the unit, before the analysis. Temperature is taken to the north of 20°N. Note that changing the domain to the north of 45°N does not vary the essential result (not shown).
Figure 9a shows the variance explained by the first 10 modes. The first mode explains 68% of the variance, and heterogeneous maps of the first mode for ablation rate and summer temperature are shown in Figs. 9c and 9d, respectively. While the second mode exhibits interesting dipole patterns over the GrIS for both variables, the explained variance is 14%, and the physical significance is not obvious from the spatial pattern (not shown). Here we focus only on the first mode. The expansion coefficients of the first mode are shown in Fig. 9b, which loosely follow the order listed in Table 4. By regressing the expansion coefficients onto the temperature and sea ice, we obtain the corresponding spatial patterns for the twenty-first-century change and late twentieth-century conditions (Fig. 10). Corresponding figures between the SVD and composite analyses resemble each other (Figs. 10a and 6a, 10b and 6b, 10c and 8c). Therefore, the robustness of the results from the composite analysis in sections 3c–e is verified by multivariate analyses.
To investigate the effect of model bias on the projected normalized ablation rate change, the spatial regressions between model biases of sea ice concentration in the late twentieth century (1982–99) and Figs. 10c and 10d are computed. The observational data used in calculating the model bias are the NOAA Optimum Interpolation (OI) SST version 2 (Reynolds et al. 2002) and Met Office Hadley Centre Sea Ice and Sea Surface Temperature version 1 (HadISST1) (Rayner et al. 2003). The spatial re gression is applied only in grids that are statistically significant at the 5% level, that is, hatched areas in Figs. 10c and 10d. The relation between the expansion coefficient of the first SVD mode and the spatial regression coefficient is shown in Fig. 11a. The correlation for the present-day winter bias (0.86) is stronger than that for the summer bias (0.67). The importance of this bias in winter sea ice distribution over that in summer is not fully understood but is possibly related to the fact that sea ice insulates heat flux from the ocean to the atmosphere more effectively during winter (Manabe and Stouffer 1979). Figure 11a indicates that most models including the ensemble mean (zero in abscissa) underestimate sea ice concentration in the hatched area of Fig. 10c, which likely results in underestimation of the ablation rate increase that is physically associated with sea ice feedbacks. Indeed, the ensemble mean underestimates the sea ice concentration in a large area of the Arctic Ocean, while it overestimates the southern ice edge (Figs. 11b–d). The underestimated sea ice concentration in the late twentieth century along the eastern coast of Greenland and in Baffin Bay (mostly Baffin Bay for winter) suppresses potential warming associated with sea ice feedbacks in GCMs. The large areal impact of a small opening of sea ice cover on air temperature was also reported from observations (Serreze et al. 2009). A similar approach can be taken for the temperature pattern in Fig. 10b. We examined model biases with respect to the National Centers for Environmental Prediction/Department of Energy Global Reanalysis 2 (NCEP-2) (Kanamitsu et al. 2002) and the 40-yr European Centre for Medium-Range Weather Forecasts Re-Analysis (ERA-40) (Uppala et al. 2005), and the results show that the models spread in both directions of positive and negative biases (not shown). The analysis conducted here is probably of more value as a guide to the significance of reducing model bias in particular regions rather than as a predictive tool. The suggestion made here through the statistical analysis, however, must be verified by sensitivity experiments with models. These results strongly suggest that accurate simulation of the present-day climate and understanding the influence of relevant processes on a range of projections are crucial.
While previous studies have shown only the spread of the future GrIS SMB projections, little attempt has been made to touch upon the source of the spread. In the current study, processes behind the spread are systematically investigated across the models. The method adopted here avoids the primary effect of coarse resolution in GCMs, and physically credible mechanisms behind the spread are proposed. Nevertheless, several potentially important processes for SMB change are not taken into account in the current study.
Unlike energy-balance-based SMB models, the effect of surface meltwater on albedo (Bougamont et al. 2005; Robinson et al. 2010; Tedesco et al. 2011) is not explicitly represented in our SMB calculations. The incorporation of this effect requires energy flux fields at much higher resolution than currently available in GCMs or online computation of the SMB on GCM subgrids. Therefore, it is not practical to apply this effect to the multi-GCM analysis. Nevertheless, the temperature index (TI) method used in the current study is constructed from observations and thus implicitly includes such an effect, although uncertainty remains under the case of much greater warming in the future.
The elevation of the ice sheet would be lowered as melting proceeds, which would then affect the temperature on the GrIS through the lapse rate. Surface albedo would also change as the ice sheet shrinks, which would also affect the temperature on the GrIS. The impact of these elevation–temperature and ice albedo feedbacks was investigated by Ridley et al. (2005) and Vizcáino et al. (2008). Ridley et al. shows that the elimination of these feedbacks does not affect the Greenland average temperature until about 100 years into the 4 × CO2 experiment, although significant local impact is reported after 300 years of the integration. Vizcáino et al. similarly shows that these feedbacks become important only after about 400 years in the 3 × CO2 stabilization experiment. While both feedbacks can be expected to be of importance on longer time scales, the estimation and understanding of SMB changes for the coming centuries is still useful even under the restriction of a fixed ice sheet topography. Ice thickness changes have to be large enough for the elevation–temperature feedback to be important, and the ice albedo feedback only becomes important when large amounts of the former ice sheet area are transformed into tundra, which takes a long time. Greve et al. (2011), on the other hand, showed that the estimated future GrIS mass loss varies significantly between the two ice sheet models (roughly 20 and 60 cm equivalent sea level rise after 300 years under the A1B ensemble mean forcing), and attributed it to different parameterizations for basal sliding and surface melting. One possible reason for the difference from our estimate might also be errors in present-day GrIS topography simulated by the ice sheet models. These uncertainties need to be explored in the future.
The effect of meltwater discharge from the ice sheet on ocean circulation is neither included in most GCMs nor in our SMB calculations. It is difficult to assess the significance of this effect because there is a spread between studies in the magnitude of AMOC weakening induced by meltwater discharge (Ridley et al. 2005; Swingedouw et al. 2007; Hu et al. 2009; Vizcáino et al. 2010; Goelzer et al. 2011).
As stated in section 2, many studies used the positive degree-day (PDD) method instead of the TI method in calculating ablation. The PDD takes into account that higher temperatures are accompanied by a longer period of melting and gives greater ablation than TI at higher temperatures. Additional calculations using the PDD scheme, with parameters as in Abe-Ouchi et al. (2007), confirm such a behavior (Fig. 12). Therefore, the PDD widens the spread of the projection. However, the PDD tends to show a much larger sensitivity to climate change than a more physically based energy balance estimate (van de Wal 1996; van der Veen 2002; Bougamont et al. 2007; Robinson et al. 2010), and the superiority of the PDD over the TI method has not been demonstrated. The uncertainty in these schemes remains to be investigated.
The GrIS SMB change for the coming centuries is estimated using output from 25 climate models, and the sources of spread in the projections are investigated. Three scenarios—SRES A1B, A2, and B1—are considered with a particular focus on the A1B scenario. In all scenarios, ablation dominates over snowfall changes. In the A1B scenario, about 60% of the intermodel spread in the ablation rate change during the twenty-first century is explained by the global annual-mean temperature change. While causality is not demonstrated in this multi-GCM analysis, other factors relevant to the spread are uncovered based on composite and correlation analyses. These are the bias in the GrIS summer surface air temperature in the present-day simulation, and the degree of AMOC weakening and sea ice reduction near Greenland during the twenty-first century. The reduction of sea ice is closely linked to the simulated late twentieth-century sea ice distribution. While the quantification of these effects requires coordinated multi-GCM sensitivity experiments, the proposed processes are physically plausible and have not been investigated systematically across the models previously. In addition, the insight gained from this study would assist in the interpretation of model differences in future projections of ice sheet mass balance, including ice dynamics.
The impact of model biases on the ensemble projection is evaluated by making use of the relation between simulated late twentieth-century climate conditions and the future SMB change. The analysis suggests that the underestimation of sea ice concentration near Greenland—the east coast and Baffin Bay in particular—likely contribute to the underestimation of the future ablation rate increase that is associated with sea ice feedbacks. While this statistically derived suggestion needs to be reexamined in sensitivity experiments with models, such an analysis provides useful information for future model development and interpretation of ensemble projections.
An important implication of the current study is that accurate knowledge of the fate of the Greenland Ice Sheet requires a comprehensive understanding of the climate system, including processes determining global mean temperature response (i.e., clouds and ocean mixing), sea ice, and ocean circulation. In addition, more effort is needed to reduce the bias in the present-day high-latitude simulations. In the meantime, relations between future SMB and present-day climate simulation, such as that established in the current study, would help narrow the range of uncertainty in the future sea level estimate, especially with greater length and variety of observational data. Moreover, such a link should be explored not only for the climatological seasonal average but also in the observed annual cycle and trends as done for other variables (e.g., Tsushima et al. 2005; Knutti et al. 2006; Hall and Qu 2006; Boé et al. 2009). Furthermore, the coordinated sensitivity experiments would be useful for verifying and quantifying the impact of individual physical processes on the ice sheet surface mass balance.
We acknowledge the modeling groups, the Program for Climate Model Diagnosis and Intercomparison (PCMDI) and the WCRP’s Working Group on Coupled Modelling (WGCM), for their roles in making available the WCRP CMIP3 multimodel dataset. Support of this dataset is provided by the Office of Science, U.S. Department of Energy. We also acknowledge NCEP2 and NOAA OI SST version 2 data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, from their Web site (http://www.esrl.noaa.gov/psd/); ERA-40 data provided by ECMWF; and HadISST1 data provided by the Met Office Hadley Centre. Technical aspects of the SVD analysis were introduced by H. Shiogama, and assistance for the PDD calculation was provided by F. Saito. Discussions with T. Sueyoshi, and comments by J. Gregory and three anonymous reviewers are greatly appreciated. W.-L. Chan kindly proofread the manuscript. We would like to express thanks for access to the computational resources of the supercomputer system (HITACHI SR-11000) at the University of Tokyo. The developers of the freely available software, Ferret and NCL, are also acknowledged. This work was supported by the Global Environment Research Fund (S-5) of the Japanese Ministry of the Environment.