We study the implications of a recently published mass balance of the Greenland ice sheet (GrIS), derived from repeated surface elevation measurements from NASA′s ice cloud and land elevation satellite (ICESat) for the time period between 2003 and 2008. To characterize the effects of this new, high-resolution GrIS mass balance, we study the time-variations of various geophysical quantities in response to the current mass loss. They include vertical uplift and subsidence, geoid height variations, global patterns of sea level change (or fingerprints), and regional sea level variations along the coasts of Greenland. Long-wavelength uplifts and gravity variations in response to current or past ice thickness variations are obtained solving the sea level equation, which accounts for both the elastic and the viscoelastic components of deformation. To capture the short-wavelength components of vertical uplift in response to current ice mass loss, which is not resolved by satellite gravity observations, we have specifically developed a high-resolution regional elastic rebound (ER) model. The elastic component of vertical uplift is combined with estimates of the viscoelastic displacement fields associated with the process of glacial-isostatic adjustment (GIA), according to a set of published ice chronologies and associated mantle rheological profiles. We compare the sensitivity of global positioning system (GPS) observations along the coasts of Greenland to the ongoing ER and GIA. In notable contrast with past reports, we show that vertical velocities obtained by GPS data from five stations with sufficiently long records and from one tide gauge at the GrIS margins can be reconciled with model predictions based on the ICE-5G deglaciation model and the ER associated with the new ICESat-derived mass balance.
The Greenland ice sheet (GrIS) stores an ice mass equivalent to a sea level rise of ∼7.3 m (Lemke et al. 2007) and therefore constitutes one of the major reservoirs of fresh water on Earth. In the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, it was concluded that a major limitation in predicting future sea level changes is the contribution from the continental ice sheets (Solomon et al. 2007). It is therefore of utmost importance to constrain the present-day mass balance of the GrIS, and hence its potential contributions to the future global sea level rise. By determining the GrIS role in the global climate system, more reliable future scenarios can be obtained. A review of the published estimates of the mass balance of the GrIS (and of the associated uncertainties) since the 1960s has been recently presented by Alley et al. (2010, see their fig. 2). The time intervals covered by the published mass balance estimates varies from few decades to less than one year, in some cases. The mass balances appear to broadly follow the time-evolution of the Greenland coastal temperature and, since the nineties, they are systematically negative. The maximum amplitude is -’267 ± 38 Gt yr-’1 for the year 2007 (Rignot et al. 2008b) (a loss of 360 Gt yr-’1 is equivalent to an average sea level rise of 1.0 mm yr-’1). As reviewed by Alley et al. (2010) (see references therein), the methods employed to constrain the mass balance of the GrIS have so far included pure assessment, surface mass balance techniques, laser and radar altimetry and, recently, gravity observations from the NASA/DLR Gravity Recovery and Climate Experiment (GRACE) satellites (see the overview of Tapley et al. 2009).
Despite the relatively low horizontal spatial resolution (300– 500 km, see Slobbe et al. 2009), GRACE provides the possibility of obtaining monthly mass time-series for the harmonic coefficients of the gravity field (Wouters et al. 2008) and to detect possible accelerations of the GrIS mass loss on relatively short observation periods (less than of one decade, see Velicogna 2009; Schrama & Wouters 2011). Published estimates of GRACE-only solutions for the average GrIS mass balance during the last decade, listed by Schrama et al. (2011) vary between -’177 ± 6 Gt yr-’1 for the period 2003–2009 (Luthcke et al. 2010) and -’230 ± 33 Gt yr-’1 for the period 2002–2009 (Velicogna, 2009). A significant part of the error on the GRACE-only mass loss estimates stems from the glacial-isostatic adjustment (GIA) correction, which contributes a secular gravity signal to the GRACE observations. Efforts have been made to control the amplitude of this correction by forward GIA modelling, using different scenarios for the viscosity profile of the mantle and deglaciation history of the late-Pleistocene continental ice-sheets (see e.g. Velicogna & Wahr 2006; Wouters et al. 2008). Recently, there have been efforts to combine GRACE and GPS data of both global (Wu et al. 2010) and regional (Ivins et al. 2011) nature when faced with ice loss from ice caps and sheets.
In a recent study, Sørensen et al. (2011) have presented a new mass balance estimate for the GrIS during the period (2003–2008), derived from repeated surface elevation measurements from NASA′s Ice Cloud and land Elevation Satellite (ICESat, see http://icesat.gsfc.nasa.gov/). ICESat is a laser altimetry mission which operated from 2002 to 2009, providing a unique data set for cryosphere studies. For ICESat, the discontinuous temporal coverage (2–3 measurements campaigns of ∼35 d yr-’1) is compensated by a substantial gain of in spatial resolution (∼5 km) with respect to GRACE and by a reduced sensitivity to the amplitude of the GIA correction (Gunter et al. 2009). Furthermore, the altimetry alone cannot provide an estimate of the mass balance (Sørensen et al. 2011), since for all altimeter systems the conversion of volume to mass changes requires modelling of the firn dynamics and of the surface ice densities (Schrama et al. 2011). Depending on the method employed to estimate the elevation change from the quasi-repeat paths and crossover points of altimetry, for the period 2003–2008 the ICESat mass balance proposed by Sørensen et al. (2011, see their table 2) varies from -’191 ± 23 (model M2) to -’240 ± 28 Gt yr-’1 (model M3). These estimates are in agreement with previous results obtained for the GrIS using other remote-sensing methods during the same time span (see, in particular, the recent summary of GRACE results by Schrama et al. 2011) and support the merits of recent ICESat-GRACE comparison studies relative to Antarctica (Gunter et al. 2009).
The spatial pattern of the M3 mass balance proposed by Sørensen et al. (2011), mapped on a quasi-regular geodesic grid, is portrayed in Fig. 1(a). In this work, the high-resolution surface load associated with this melting scenario will be employed for providing estimates of time variations of various geophysical quantities. These estimates are essentially based on the solution of the sea level equation (Farrell & Clark 1976) and on the application of well-established Green′s function methods (Farrell 1972) that permit to capture the shorter wavelength features of the regional uplift patterns (Khan et al. 2007). Predictions include long- and short-wavelength vertical deformations in response to unloading, regional sea level variations along the coasts of Greenland, and global sea level patterns (or fingerprints) associated with the current mass loss of the GrIS. In the last decade, significant efforts have been devoted to the observation and interpretation of GPS signals possibly associated to the present-day vertical movements across Greenland (Wahr et al. 2001; Khan et al. 2007, 2008, 2010a; King et al. 2010), with the aim of weighing the relative importance of GIA and elastic rebound (ER) in response to current ice melting. Using the observations from the few stations that can currently provide a reliable trend of vertical uplift (Khan et al. 2007, 2008), we will address the problem of reconciling the GPS vertical uplift data with GIA and ER modelling.
The paper is organized as follows. The methods are described in Section 2. In the results Section 3, we present the global elastic deformations (Section 3.1), the pattern of global and regional sea level variations (3.2), and the regional elastic uplift of the GrIS (3.3) in response to the current mass loss derived from ICESat observations. After a discussion of the geophysical implications of the results (Section 3.4), we draw our conclusions.
Since Greenland is recognized as a tectonically stable region with low seismicity (Chung & Gao 1997), we can assume that the present-day uplift originates from two major components acting simultaneously (e.g. Khan et al. 2008). The first is the ER associated with the recent ice mass loss across the GrIS, which has been observed by geodetic techniques in the last few decades. In view of the very short timescales involved, it is generally assumed that the viscous component of this deformation is negligible (Spada et al. 2011a). The second, associated with GIA, results from the viscoelastic response of the bedrock to the evolution of the GrIS and of the surrounding ice complexes during the last glacial-interglacial period (a complete summary of state of the art GIA theory is presented in the recent report of Whitehouse 2009). These two geophysical processes and the methods employed for their numerical simulation are synthetically described in Table 1. Though ER and GIA are characterized by distinctly different timescales and spatial patterns, we know from previous analyses that their amplitudes may be locally comparable, and therefore they are both important for a correct interpretation of GPS observations in Greenland (e.g. Khan et al. 2008). In what follows, efforts will be made to estimate modelling uncertainties of the GIA and of the ER component of the total deformation. These uncertainties are related with a limited knowledge of the driving mechanism (i.e. the history of ice thickness variation) and with the response of the solid Earth (i.e. its elastic and viscoelastic structure).
The GIA component of sea level change is evaluated solving the sea level equation, introduced by Farrell & Clark (1976). It readsSpada & Stocchi 2006). Assuming a Maxwell viscoelastic rheology, eq. (1) is solved for S to a maximum harmonic degree ℓmax by the the pseudo-spectral iterative method (Mitrovica & Peltier 1991) using code SELEN (Spada & Stocchi 2007). Peltier, 2001). Vertical displacement can be written as Spada & Stocchi 2006). Apart from the notation, these definitions of N and G are matching those by Simon et al. (2010). In this study, we are mainly concerned with the trend (i.e. the time-derivative) of the above quantities, which in the body of the paper will be denoted by , , and , respectively, and assumed to be constant over the time period covered by the ICESat observations (∼2003–2008).
The ER problem for Greenland will be solved using both a global and a regional—GrIS scale—approach (see Table 1). In the following, these will be referred to as GER and RER approaches, respectively. The elastic LDCs required for these methods have been obtained from the atmospheric pressure loading service (here after abbreviated to as APLO). They pertain to a spherically symmetric, non-rotating, elastic and isotropic (SNREI) earth model with PREM structure (Dziewonski & Anderson 1981). They are expressed in the reference frame of centre of mass of the system (Earth + Load) and account for the Earth′s elastic compressibility. The LDC for vertical displacement (hAPLO), is shown in Fig. 2(a) as a function of degree ℓ. The asymptote for ℓ↦∞ is almost attained, here, for ℓ≈ 103. For a planet of radius a, mass m and gravity γ, , where λ and μ are the Lamè constants of the surficial layer of the Earth (Farrell 1972).
While GER relies upon the solution of the full sea level equation, in the RER method we simply perform a convolution between the elastic Green′s functions and each of the elementary disc-shaped elements that discretize the mass budget of the GrIS (see Fig. 1a). Apart from the different Green′s functions employed, the RER approach thus follows closely the method adopted by Khan et al. (2007, 2008) to evaluate the vertical movements at GPS sites in response to the melting of near-field GrIS glaciers. For a single ice element, vertical displacement isSpada et al. 2011b). A similar expression holds for the variation of the geoid height G(α), with hℓ substituted by 1 +kℓ. Once eq. (4) is evaluated for a unit ice thickness variation, the total displacement U(θ, λ) corresponding to a regional ensemble of disc loads is obtained by rescaling the solution and taking advantage, by means of simple geometrical methods, of the axial symmetry of eq. (4) relative to the centre of the load. Fig. 2(b) shows U(α) for an instantaneous melting episode with ΔH=-’1m. The load amplitude is β= 0.025°, matching the size of the 5 km × 5 km ice elements that compose the GrIS M3 mass balance. Displacements are shown for increasing values of ℓmax until the oscillations in the periphery of the load (see inset), have disappeared (the final solution corresponds to ℓmax=ℓ∞≡ 105). Uplift is attained across a broad area surrounding the load (approximately for α≤ 16β); for larger values of α, the load is surrounded by a broad peripheral subsidence. As a rule of thumb, a melting episode with ΔH=-’1.5 m produces, if hAPLO LDCs are employed, a vertical uplift of ≈1 mm beneath the load.
Differently from GER, in the RER method, deformation only results from the direct effect of ice melting; meltwater loading and the gravitational interaction between the GrIS and the surrounding oceans are not included in modelling. The GER and the RER approaches have a complementary character. Being ‘gravitationally self-consistent’ (Wu & Peltier 1982), the GER method allows for estimates of global and regional deformations and gravity changes, and therefore provides, in principle, reliable predictions of sea level variations at any spatial scale. However, since the computation of short-wavelength deformations by the GER method may be computationally very expensive, for regional and local predictions the RER approach is by far more convenient and, as we will discuss in Section 3 below, largely consistent with GER on the spatial scale of the GrIS.
Global elastic deformations in response to current mass loss
Figs 3(a)–(c) show the rates of vertical displacement (), of relative sea level change () and sea surface variation (), respectively, based on the GrIS mass budget M3 and obtained by the GER method to degree ℓmax= 128. The fields are mapped on a quasi-regular geodesic grid with spacing of ∼25 km, obtained by the icosahedron method of Tegmark (1996). In the range of wavelengths considered, the maximum rates of elastic uplift in Fig. 3(a) amount to ∼10 mm yr-’1 and are positively correlated with the regions of coastal ice thinning, where mass loss is reaching its largest amplitude according to the mass balance M3 shown in Fig. 1(a). However, the local vertical uplift is so strong that the whole GrIS is uplifted, even in the central part where the ICESat-based mass balance indicates ice accumulation and our predicted velocities are in the range of 2–4 mm yr-’1. Due to the limited spatial resolution of this uplift map, it is not possible to detect areas of local subsidence, which are however suggested by the localized ice thickening along the NE coastline of Greenland in Fig. 1(a). The rate of sea level change , shown in Fig. 3(b), closely (but not exactly) mirrors the map of , and shows a complex regional pattern along the coasts of Greenland. It reaches its largest amplitude (∼-’14 mm yr-’1) in the proximity of the Helheim Glacier (HG, see Fig. 1b). Fig. 3(c) shows that the sea surface, as seen from the Earth′s centre of mass, is collapsing across Greenland. Furthermore, the pattern of is clearly characterized by less short-scale variations compared with and , which is a consequence of the larger energy contained in the low-degree portion of the gravity field spectrum.
In Fig. 4, a comparison is drawn between the uplift pattern obtained by the ICESat mass balance M3 of Sørensen et al. (2011) and that inferred by observations from the GRACE satellites. Among the various GRACE solutions so far published (e.g. Chen et al. 2006; Luthcke et al. 2006; Ramillien et al. 2006; Velicogna & Wahr 2006), here we have used the one by Wouters et al. (2008), who have employed 58 monthly GRACE observations between 2003 February and 2008 January, a time span that substantially overlaps the period of the ICESat observations considered by Sørensen et al. (2011). The GRACE mass balance estimate is derived from the level-2 data, from which de-aliasing products such as atmospheric and non-tidal ocean loading are subtracted. In the GRACE data such loading has a direct mass effect while in the ICESat altimetry it will have the indirect effect of bedrock movement. Ocean and atmospheric loading is not taken into account in the ICESat derived mass balance used here, because the Earth′s response to these is a high frequency signal, and it is assumed to have a small (if any) influence on the trend in surface elevation (van Dam et al. 1997).
According to the GRACE results of Wouters et al. (2008), the GrIS has lost mass at an average rate of 179 ± 25 Gt yr-’1, an amount somewhat smaller than the M3 ICESat estimate of -’240 ± 28 Gt yr-’1 (Sørensen et al. 2011). We note that this GRACE solution agrees better with the mass budget based on method M2 by Sørensen et al. (2011), namely -’191 ± 23 Gt yr-’1. Since the GRACE mass balance is obtained by an iterative procedure in which Stokes coefficients of the gravity potential are inverted to degree ℓmax= 60, in Fig. 4(a) the ICESat-derived uplift rates have been truncated to the same maximum degree and the same sets of APLO vertical LDCs are employed for comparison. For ICESat, the truncation to degree ℓmax= 60 implies a substantial loss of information, as it can be appreciated by a comparison between Fig. 4(a) and Fig. 3(a) (obtained with ℓmax= 128). The misfit between ICESat and GRACE, portrayed in Fig. 4(c), attains its minimum in the bulk of the GrIS and in the northeast, where the differences do not exceed 1 mm yr-’1. However, the GRACE mass balance does not resolve in detail the ice loss captured by ICESat in the area of the Southeast Glaciers (SG) nor the very localized melting in the area of the Jakobshavn Isbræ (JI), in the northwest of the GrIS (see also Fig. 1). In these regions, the predicted ICESat uplift rate exceeds the GRACE solution by 1– 3 mm yr-’1et al. Considering the completely different nature of the two data sets and the different inversion methods, the results of Fig. 3 support the view that GRACE and ICESat observations are providing, in addition to comparable gravitational signatures (Sørensen 2010), consistent scenarios for the long-wavelength elastic uplift patterns during the time span 2003–2008. Differences between the mass balances obtained, and consequently discrepancies in the vertical uplift predictions as shown in Fig. 3 can be attributed to (i) undersampling of the southeastern part of the ice sheet by ICESat (Sørensen 2010), (ii) to the intrinsic limitations that all satellite altimetry systems face along the coastal margins because of the rough topography (Schrama et al. 2011), (iii) to the different temporal and spatial resolutions of GRACE and ICESat and (iv) or to the not exact coincidence of the time periods of the ICESat and GRACE time-series used for the comparison of Fig. 4.
Sea level fingerprints
The GER approach based on eq. (1) allows for a ‘gravitationally self-consistent’ evaluation of the global sea level variations associated with the M3 mass balance model of the GrIS. The sea level change pattern or ‘fingerprint’ (Mitrovica et al. 2001, 2009; Riva et al. 2010) associated with this mass balance solution is shown in Fig. 5. With a total mass budget of -’240 ± 28 Gt yr-’1 in the period 2003–2008, according to our computations the GrIS model M3 would provide, on a rigid and non-self-gravitating Earth, an eustatic sea level variation at a rate of +0.67 ± 0.08 mm yr-’1. This estimate agrees well with results found by others. For example Wouters et al. (2008) found an estimate of 0.5 ± 0.1 mm yr-’1 based on GRACE for the time period between 2003 February and 2008 and January. The estimate is also consistent with previous GRACE satellite gravity observations (Velicogna & Wahr 2006) which indicate a mean mass loss between 2002 and 2006 equivalent to a uniform sea level rise of 0.5 mm yr-’1. Although based on different observations (i.e. tide gauge records) and representative of a secular time scale, Mitrovica et al. (2001) also proposed a comparable estimate of +0.60 ± 0.15 mm yr-’1. Fig. 5 clearly illustrates significant deviations from eustasy, with maximum values that exceed the average eustatic value (white contours) by ∼15 per cent in the regions located in the far field of the GrIS (e.g. in the middle of the Pacific Ocean). This agrees with the recent results of Riva et al. (2010), who estimated a gravitationally driven sea level increase of ∼20 per cent greater than the average (eustatic) value at low-latitudes (this slight difference is likely to result from the effects of Earth rotation, which are included in Riva et al. 2010). In Fig. 5(b), the same total mass budget employed in 5(a) has been uniformly distributed over the GrIS, therefore following the same approximation of Mitrovica et al. (2001). The corresponding mass balance, henceforth referred to as M3U, causes the same eustatic sea level variations of M3. Differences observed with respect to the (normalized) sea level pattern of fig. 1(b) of Mitrovica et al. (2001) can be attributed to the absence of rotational feedbacks (Milne & Mitrovica 1998) in this elastic implementation of the sea level equation, and possibly to the different spatial resolutions employed [Mitrovica et al. (2001) use ℓmax= 512, while in Fig. 5ℓmax= 128)]. Differences between the maps in Figs 5(a) and (b) are evident (and expected) in the near field of the GrIS. Despite the significant spatial heterogeneity of the actual mass loss distribution employed in Fig. 5(a), associated with the east–west bipolar structure of the mass balance of model M3, differences between the two fingerprints in the remote far field cannot be appreciated at the 0.1 mm yr-’1 level. Similar results were observed by Bamber & Riva (2010).
The geophysical significance of the results shown in Fig. 5 can be better established by a comparison with the secular rates of sea level change attributed to GIA in key areas as the Mediterranean basin. Fig. 6 shows the details of the sea level trends across southern Europe and the Mediterranean Sea, adopting the same melting scenario as in Fig. 5(b). The values vary smoothly within the basin (approximately from 0.4 to 0.6 mm yr-’1 from NW to SE), due to the very long wavelength character of the far-field sea level variations driven by the mass loss in Greenland. Although a precise estimate is premature, considering the uncertainties in the mantle viscosity profile, according to Tsimplis et al. (2011), the GIA component of the sea level variations currently observed at tide gauges located in bulk of the Mediterranean Sea amounts to ∼0.5 mm yr-’1. This estimate, based on the global ice sheet reconstruction ICE-5G(VM2) of Peltier (2004), is consistent with independently derived models of deglaciation (see the very recent work of Lambeck et al. 2011). The GIA secular sea level trend is comparable with the values shown in Fig. 6 in the same region, which is only based on the M3 ICESat mass budget. For tide gauge sites located along the continental coasts of the Mediterranean Sea, such as Trieste (northeast Italy), the current melting of Greenland largely dominates the GIA effects, which according to Tsimplis et al. (2011) may be close to 0.1 mm yr-’1. This is consistent with the GIA modelling results by Stocchi & Spada (2009) and Lambeck et al. (2011). Since here we only consider the effects of mass loss from the GrIS during less than one decade, neglecting contributions from Antarctica (e.g. Rignot et al. 2008a) and small glaciers (e.g. Meier et al. 2007), these estimates should be taken cautiously since they may be not representative of secular sea level trends. However, they indicate that current ice melting of the GrIS is contributing a large fraction of the current sea level variations observed at tide gauges in the Mediterranean. The distinct regional character of the sea level fingerprints shown in Fig. 6 should be taken into account for future sea level predictions in this region, which are often based on the eustatic approximation (e.g. Lambeck et al. 2011).
Regional elastic deformations
The GER method that we have so far employed for modelling ER of Greenland is based on a solution of the sea level equation in which only large to intermediate wavelengths are taken into account (ℓmax= 128). However, the detailed surface load variations provided by ICESat observations, surface density and firn compaction modelling on the 5 × 5 km grid of Fig. 1(a) may have the potential of affecting vertical movements (and consequently regional sea level variations) on a much smaller spatial scale. To improve the spatial resolution of our regional predictions, we have first set up a RER experiment based on eq. (4) and mapped the rate of displacement in the same range of wavelengths considered so far, using the APLO LDCs shown in Fig. 2(a). The results, shown in Fig. 7(a), are matching well those obtained by the GER method (Fig. 3). The difference between the GER and RER results for , shown in Fig. 7(b), is positive across the whole GrIS and reaches its maximum off the coast segments which are more subject to rapid mass loss (and where consequently vertical uplift is maximum). This extra regional uplift in the range of 0.1–0.3 mm yr-’1 originates from the sea level fall associated with the decreased gravitational attraction between the melting ice masses and the oceans surrounding the GrIS. This effect, associated with the self-gravitation of the ocean masses, is accounted for in the sea level equation but totally neglected within the RER modelling. Since its amplitude is substantially negligible compared to the values predicted along the GrIS margins, the results of Fig. 7 indicate that the RER method provides, on these regional scales, an accurate estimate of the uplift pattern despite the significant simplifying assumptions in the formulation (see Section 2). This finding has important practical implications, since, at the same ℓmax, the RER approach is computationally more convenient than the physically self-consistent (but expensive) GER method. As shown recently for Antarctica by Simon et al. (2010), the ocean loading effect on the vertical crustal motion is larger when viscoelasticity is taken into account.
The short-wavelength features of the elastic component of the vertical velocity field and of the rate of geoid variation are shown in Fig. 8 on the GrIS spatial scale. Detailed views are shown in the snapshots of Fig. 9. In both figures, we show results of RER computations to harmonic degree ℓ∞= 105 based on the hAPLO LDCs (the corresponding response function is shown in Fig. 2b). We note that no broad-scale subsidence is observed, even in the central portion of the GrIS where our ICESat mass balance indicates mass increase during the period 2003–2008 (see Fig. 1). This is due to the broad range of the elastic uplift that, according to the results of Fig. 2, is surrounding the centres of mass loss in spite of their small spatial wavelength (a relevant discussion on this character of the elastic response to surface loads is given by Stocchi et al. 2005). This long-range uplift prevails on the subsidence induced by ice accumulation in the central portion of the GrIS. According to Fig. 8(a), in the bulk of the GrIS the elastic uplift varies in the range between 2 and 4 mm yr-’1, comparable to that estimated by the low-resolution GER method in Fig. 3(a). Similarly, the large-scale rate of geoid subsidence is comparable to that of the sea surface considered in Fig. 3(c). However, these highly detailed maps allow resolution of quite realistically isolated features of the vertical velocity field, whose amplitudes significantly exceed that obtained by the low-resolution GER analysis. Maximum predicted rates are mm yr-’1 and mm yr-’1, observed in the vicinity of the Kangerdlugssuaq Glacier (KG). These improved estimates may have an impact on the interpretation of GPS data from stations located close to glaciers associated with rapid mass loss. This is the case, for example, of the GPS site of KULU, located in the SG region (see also Fig. 9e). In fact, as shown in Table 3, our RER computations at KULU deviate quite significantly (by about ∼2 mm yr-’1) from those obtained by the GER method. At the other GPS sites relevant for this study, located at larger distances from the main drainage basins, the differences are less significant, to indicate that these sites are mostly responding the large scale uplift of the GrIS, which results from the superimposition of the contributions of all elementary ice elements of Fig. 1(a) and cannot be directly attributed to individual localized glaciers. Due to the different mass balances employed, a direct comparison of our RER computations with the elastic computations of Khan et al. (2007, 2008) is not straightforward (but it will be attempted in Section 3.4). However, since in Figs 8 and 9 we account for the whole mass balance, we expect that ICESat-derived vertical rates of displacement may exceed those of Khan et al., who have only used the mass balance from the major drainage systems.
Due to the increased spatial resolution of our RER computations, localized areas of subsidence appear very clearly in Fig. 8 on the scale of a few tens of kilometre, particularly along the NE coasts of Greenland. The existence of these subsiding areas could be deduced from Fig. 1(a), which reveals localized spots of positive mass balance (Sørensen et al. 2011). In particular, subsidence occurs in the area of the outlet glacier of Storstrømmen (see detail in Fig. 9b), which is recovering from a previous surge that lasted until 1978 (Bøggild et al. 1994; Reeh et al. 1994) and where according to our computations the rate of subsidence is reaching an amplitude of ≈2 mm yr-’1. However, due to the resolution of the ICESat data, the complete mass loss signal might not be captured in the outlet glaciers, meaning that the elastic uplift signal might be higher (especially near the glacier front) than what can be seen from ICESat. Another area of subsidence is predicted in the vicinity of the Flade Isblink glacier, the largest ice cap in Greenland (Fig. 9a). As discussed by Sørensen (2010), at the Flade Isblink ice cap, solution M3 suggests a complex pattern of elevation change. While the northwestern part of the ice cap has thickened with maximum values of more than 1 m yr-’1, the southeastern part shows a thinning of up to 0.5 m yr-’1. This supports recent results of Kelly & Lowell (2009), who indicate that at present Flade Isblink is growing, though further research is under way by Center for Remote Sensing of Ice Sheets (CReSIS) (John Paden, personal communication 2011) to examine this growth. According to our RER computations, the subsidence rate in the Flade Isblink can be estimated at the level of ∼0.5 mm yr-’1. Note that these estimates of subsidence in the NE of the GrIS only account for the current elastic uplift. In the same regions, strong regional scale uplifts associated with the GIA component of vertical displacement would overwhelm these areas of subsidence, at least according to our computations based on model ICE–5G of Peltier (2004) (see Fig. 10a).
The modelling results presented in the previous sections indicate that the elastic uplift of Greenland in response to current mass loss is characterized by an heterogeneous pattern (see, in particular, Fig. 8). Localized, glacier-scale vertical movements at rates as large as tens of millimetres per year are superimposed on a background GrIS-scale uplift of a few millimetres per year that ultimately results from the interference of the localized sources. Since Greenland is considered a stable continental region with low seismicity (Chung & Gao 1997), GIA represents the second major cause of vertical deformation in the area (seismicity is mainly localized along the coastlines and possibly associated with the GIA process itself, see Chung 2002). To estimate upper and lower bounds of the response of Greenland to the northern hemisphere ice sheets fluctuations (i.e. the GIA component of vertical deformation), we have computed, via eq. (1), the rates of vertical displacement expected across the GrIS according to a set of various plausible scenarios of the global melting history since the Last Glacial Maximum and the Earth rheological profile. This approach, similar to that employed by Velicogna & Wahr (2006) and Simpson et al. (2011) differs from previous studies in which a single model of deglaciation has been considered (see e.g. Khan et al. 2008) and can tentatively provide the range of uncertainly associated with GIA modelling. These uncertainties arise as a consequence of the different observational constraints used to reconstruct the history of melting, different assumptions about the rheological profile of the mantle and different numerical implementations of the sea level equation. We also remark that we did not attempt any simultaneous inversion of GIA parameters based on classical relative sea level (RSL) observations and space geodetic methods, as done for example by Paulson et al. (2007). Rather, in our GIA simulations, we have used fully a priori ice sheets scenarios and viscosity profiles, being aware that a significant modification of the viscosity profile would amount to generate a new ice history since the Last Glacial Maximum in order to match the RSL observations (Paulson et al. 2007).
Our GIA computations, obtained by an improved version of the open source code SELEN (Spada & Stocchi 2007), are shown in Fig. 10. Here we account for the sea level variations associated with the rotational feedback (Milne & Mitrovica 1998) and, since we aim to estimate the uncertainties related with the ice sheets history and mantle viscosity, we do not account for the horizontal migration of shorelines. In all computations, the deglaciation scenarios for North America, Eurasia, Antarctica and minor ice sheets are also included. As in our previous RER and GER computations, rates of displacement are computed in the reference frame with origin in the centre of mass of the Earth system (including surface loads). In Fig. 10, we show results obtained by implementing four different scenarios in our numerical code SELEN. In Fig. 10(a) we employ the ice sheets chronology ICE–5G (Peltier 2004). More specifically, we have implemented the ‘ICE–5G (VM2 L90) model version 1.2’ (1° resolution), available from the home page of Prof. W. R. Peltier. In Figs 10(b)–(d) we consider three variants (corresponding to different viscosity profiles, see below) of the ice model progressively developed at the Research School of Earth Sciences (RSES) of the National Australian University by Kurt Lambeck and coworkers (see Fleming & Lambeck 2004, and references therein). Since this model, kindly provided by K. Lambeck, is valid as of 2005, we will refer to is as to ANU05 in the following. Its Greenland component is referred to as GREEN1—first order Greenland ice model—by Fleming & Lambeck (2004). For details regarding our numerical implementation of ICE–5G and ANU05, the reader is referred to Table 1. We remark that details of our numerical implementation of the sea level equation may differ from those in the original works where these two ice models have been presented.
Models ICE–5G and ANU05 differ in several aspects, including the mantle viscosity profile, the ice volumes distribution at the Last Glacial Maximum and the time history of equivalent sea level, which represents the amount of ice water exchanged between the oceans and continental ice sheets over time. These differences are ultimately the consequence of the different sets of global relative sea level data and modern geodetic observations used by the authors to constrain the model parameters. Here, the ice histories of ICE–5G and ANU05 are implemented in SELEN (Spada & Stocchi 2007). For ICE–5G, which incorporates the Greenland ice sheet reconstruction by Tarasov & Peltier (2002), we employ an incompressible earth model characterized by a two-layer viscosity profile that approximates the multilayered profile VM2 (details are given in Peltier 2004) while for ANU05 we adopt the effective values for the lower, nominal and upper solutions of Fleming & Lambeck (2004) (viscosity values are shown in Table 2 for the four models employed here). Lower and upper viscosity solutions are not prescribed, and simply represent limits on the range of viscosities. We note that independent computations based on ICE–5G (Khan et al. 2008, see their fig. 9), show the same broad pattern of Fig. 10(a), but with significantly larger local values (in the range of several millimetres per year) and a markedly enhanced subsidence in the Baffin Bay and in the Davis Strait. As discussed below, these differences have a significant impact on the interpretation of the GPS observations. Though it was not possible to unequivocally identify the origin of these differences, it is likely that they result from a combination of effects, including the time-discretization employed to compute the time-derivatives of vertical displacements and geoid variations, the different spatial resolution (i.e. the ℓmax value) employed in modelling (L. Tarasov, personal communication 2011) and our choice to use a coarse two-layers depth-average version of the VM2 viscosity profile. We also note that predictions by Khan et al. (2008) were directly obtained from the page of SBL (the Special Bureau of Loading of the International Earth Rotation Service, see see http://www.sbl.statkart.no), devoted to the outputs of model ‘ICE–5G V1.2’. The temporal discretization of this model, however, apparently differs from that characterizing the ‘ICE–5G (VM2 L90) model version 1.2’ implemented in our numerical simulations and nowadays available from the web page of Prof. W. R. Peltier. In fact, the former includes an ice loading phase starting 122.0 kyr ago and has very recent ice thickness variations (with a significant ice thickness variation over Greenland during the last 100 yr), while the latter only describes the melting history between the Last Glacial Maximum (21 ka BP) and present time in time steps of 0.5 and 1.0 kyr, with the most recent ice thickness variation occurring during the last 500 yr. Though we did not perform a detailed analysis of the consequences of these differences on present-day vertical movements in Greenland, they could, at least partly, explain the discrepancies between the uplift pattern shown by Khan et al. (2008) and our numerical results in Fig. 10(a). These issues are going to be addressed within a sea level equation benchmark program along the lines of a recent initiative of the GIA community (Spada et al. 2011b). To facilitate the intercomparison and the validation of our results, the GIA numerical codes and input data used to obtain Fig. 10 are available from the author.
The maps in Fig. 10 show a significant variance that supports previous concerns about the effective accuracy of GIA computations over Greenland (Khan et al. 2007, 2008). For all solutions shown, local vertical velocity values across the GrIS are of the order of a few millimetres per year, somewhat less than the background (regionally averaged) RER elastic velocities shown in Fig. 9. Correcting the ICESat elevation change data for these GIA effects implies a mass budget correction of approximately 1 Gt yr-’1 (Sørensen et al. 2011), a small amount if compared to the cumulative error associated with the ICESat data inversion and firn compaction corrections (±28 Gt yr-’1 for solution M3). As it is apparent from Fig. 10, a major difference between ICE–5G and ANU05 results is the uplift pattern in the southwest of Greenland, where the ICE-5G predictions show a marked subsidence in the surroundings of the GPS site of KELY, approximately between latitudes 62° and 72°N. As pointed out by Khan et al. (2008), this subsidence is caused by the re-advance of the ice margin in west Greenland to its current location during the last few kilo-years (the subsidence at KELY was first explained by Wahr et al. 2001, by a re-advance model). This feature appears not to be present in the ANU05 model, though the effects of more recent re-advances (‘neoglacial’ periods) have been considered by Fleming & Lambeck (2004). Other significant differences between the predictions based on ICE–5G and ANU05 visible in Fig. 10 can be attributed to the characteristics of the rheological profiles employed (see Table 2), to the time history of the ice sheets surrounding the GrIS (especially Laurentia) or to the methods followed in the GrIS reconstruction. While Fleming & Lambeck (2004) have only employed constraints from Holocene RSL histories from sites along the coastlines of Greenland, Tarasov & Peltier (2002) have also used 3-D dynamic ice sheets models (e.g. Huybrechts, 1996). As discussed below (and as pointed by Khan et al. 2010a), observations of vertical movements at the GPS sites shown in Fig. 10(a) in conjunction with ER modelling can help to discriminate between competing GIA models. However, a full understanding and interpretation of some features of the uplift patterns in Fig. 10, such as those apparent in the northeastern portion of the GrIS will only be possible when sufficiently long series of data will become available from the GNET sites shown in Fig. 1(b).
In Fig. 11 we compare various model predictions with GPS observations at the GPS sites of KELY (a), THUL (b), SCOR (c), KULU (d) and QAQ1 (e), whose locations are reported in Fig. 1 (the numerical data employed to draw this figure are summarized in Table 3). Black diamonds reproduce the results of Khan et al. (2008) in their table 4 and only account for ER. Similarly, black triangles are results from our RER modelling based on the mass balance M3 shown in Fig. 1(a) and on the response function of Fig. 2(b) with ℓmax=ℓ∞. By further computations, we have verified that using model M2 of Sørensen et al. (2011), characterized by a mass balance reduction of ∼20 per cent relative to M3, would shift the RER results by ∼1 mm yr-’1 downward without changing the main conclusions of this work. The combination of the RER predictions with the four GIA scenarios discussed in Fig. 10 are displayed by an open circle (ICE–5G) and squares (the three variants of ANU05). Observed GPS rates of vertical uplift and their uncertainties (±1σ), shown by shaded rectangles, have been provided by Shfaqat Abbas Khan (personal communication 2011). Since the vertical rates caused by elastic deformation vary with time (Khan et al. 2010a) and the ICESat-derived uplift rates are based on observations between 2003 and 2008, here we use GPS rates of vertical displacement obtained using 2003–2008 GPS data that update previous results by Khan et al. (2008). They have been processed using the GIPSY OASIS 6.1 software package (Zumberge et al. 1997) and using IGS orbits, Earth orientation parameters, and clock products as described by Khan et al. (2010b). The solutions are aligned with the IGS05 frame Altamini et al. (2007). For NUUK (a tide gauge station), the rate of vertical uplift has been obtained by Khan et al. (2008) with the aid of altimetry data. Following Khan et al. (2008), we make the hypothesis that the secular trends at the sites considered in Fig. 11 are associated with the combined effects of the elastic uplift caused by the present-day mass loss and of the (steady-state) viscoelastic uplifts associated with past changes. However, short-term effects associated with crustal heterogeneity (Ivins & Sammis 1996) or shallow upper-mantle transient components of deformation (Spada et al. 2011a) cannot be discounted, even on a decade timescale. Furthermore, inter-annual effects caused by accelerations of the mass loss of the main outlet glaciers of the GrIS (Velicogna & Wahr 2006; Khan et al. 2010a) are not considered here, since the ICESat mass budget employed in this study represents a time average over the whole observation period (Sørensen et al. 2011).
Some trends in Fig. 11 are worthy of discussion. First, ER (black filled symbols) effectively appears to constitute the major component of total vertical deformation at most of the GPS sites considered. However, there are remarkable exceptions which will be discussed below. Second, the RER solution (see triangles) always exceeds (in some case very significantly, as for THUL), the Khan et al. (2008) results (diamonds). We expect that the dominance of the RER uplift rates mostly reflects differences in the mass budgets employed. To test this hypothesis, we have adopted an approach similar to Khan et al. (2008), who have employed the regional mass balances published by Rignot & Kanagaratnam (2006), based on radar interferometry observations. In particular, we have performed further RER computations in which we have retained, in the ICESat mass balance M3 of Fig. 1, only the drainage basins numbered by 10, 11, 13 and 20 in the work of Rignot & Kanagaratnam (2006), which are encompassing the glaciers KG, HG, SG and JI, respectively. The results of this consistency test, reported in Table 3, are found to be in good agreement with those obtained by Khan et al. (2008) (the only significant disagreement is possibly relative to KULU). This finding further confirms the consistency of the ICESat-derived mass balance M3 with independent observations, and explicitly indicates that the contribution of relatively small glaciers located in the vicinity of the GPS sites (and naturally included in our RER modelling) is important and in specific cases it can dominate the one associated with larger but relatively distant sources. In this respect, a particularly interesting example is that of THUL, where the RER vertical uplift rate (Fig. 11b, ∼6 mm yr-’1), largely exceeds the small value computed by Khan et al. (2008) (0.3 mm yr-’1), who have only included the major glaciers in their computations. It is likely that the vertical deformations at the GPS site of THUL are mostly associated with the mass loss at the nearby Gades Glacier and at those further to the south, where according to our modelling, uplift rates are locally as large as ∼15 mm yr-’1 (see Fig. 9f).
According to our computations in Fig. 11, with the exception of KULU predictions based on the combination RER+ICE–5G (circles) are consistently within the shaded (±1σ) region or are marginally outside. Considering the modelling uncertainties and the simplifying assumptions made, we find that this level of agreement between GPS observations and model predictions is largely acceptable. In QAQ1, the GIA process does not alter the already satisfactory fit between our RER solutions and the GPS observations while, remarkably, the GIA contribution helps to reconcile the RER solution with the GPS trend at the sites of THUL, NUUK, KELY and SCOR. This result is clearly at variance with one of the main conclusions of Khan et al. (2008), who attributed the significant discrepancies (as large as ∼3.7 mm yr-’1) between observed and predicted rates at THUL and SCOR to limitations of the GIA model (namely, ICE–5G). The misfit between all model predictions and the GPS observation of KULU is apparent. A possible explanation is that ICESat underestimates the marginal thinning occurring in the area of the Helheim glacier. ASTER and ATM data indicate a thinning up to 200 m near the glacier front during 2003–2006 (Howat et al. 2008), where ICESat tracks are not close enough to the glacier front and cannot capture the maxima of the signal. The problem is the same for several other SE glaciers. As a further remark, we note that the model predictions including the GIA component based on model ANU05 (squares) show a considerable scatter and, globally, a reduced performance compared to those based on ICE-5G. The misfit can be tentatively attributed to the inhomogeneity of the spatial distribution of the RSL observations that are used to constrain the spatiotemporal features of ANU05. For example, the lack of RSL observations from the southeast of Greenland (Fleming & Lambeck 2004) can be invoked to motivate the relatively poor fit at KULU, but this could hardly justify the misfit at KELY and NUUK, seen the relatively large number of observations that Fleming and Lambeck have employed from the coasts facing the Baffin Bay. Although it is not our purpose to perform a rigorous statistical analysis of the results in Fig. 11, the misfit diagram shown in Fig. 12 can help to classify the performances of the various models employed. The misfit curve attains its largest values when only the effect of GIA is considered (right portion of the figure) and shows a more complex behaviour when ER or ER+GIA models are considered (see left and central portion, respectively). However, in agreement with the qualitative scenario that we have drawn by Fig. 11, the best performing combination of models is RER+(ICE–5G) (M≃ 5). For the dashed curve, obtained excluding the site of KULU from the computation of misfit, for the best-fitting model M reduces to ≃2) (as a rule of thumb, a misfit M≈ 1 would indicate a good match between data and model predictions).
In Fig. 13, we show the rates of vertical uplift expected at all the GPS (GNET) sites of Fig. 1(b), according to different combinations of GIA and RER computations. Dotted curves only account for the elastic deformation, while the solid ones incorporate the GIA (ICE–5G) effects. Due to the short record period, vertical uplifts at these GPS sites cannot still be established with an acceptable level of precision. Hence, observed rates are not shown except those at the main GPS sites of Table 3. Assuming that the rate of ice loss in Greenland will remain constant in the future and that transient rheological effects (Ivins & Sammis 1996) can be neglected, the results shown in Fig. 13 can provide an estimate of the rates of vertical uplift expected during the next century (Spada et al. 2011a). According to the results of Figs 3(a)–(c) to a first approximation regional sea level variations along the coasts surrounding the GPS sites will be mainly controlled by vertical displacements, with an amplitude ≈-’dU/dt. Fig. 13 shows the relative importance of GIA compared to ER. GPS sites in the range marked by the two grey shaded rectangles in frame (a) are virtually only affected by ER. It clearly appears that among the GPS sites that currently dispose of a relatively long record period, only THUL and KELY are sensitive to the GIA component of vertical deformation.
We have investigated some geophysical consequences of the recently published estimate of the annual mass loss of the GrIS (Sørensen et al. 2011), derived from surface elevations observed by ICESat observations, firn compaction and surface density modelling. In particular, assuming a radially stratified structure for the mantle, we have provided new estimates for the crustal uplift rates associated with elastic and viscoelastic rebound, regional and global sea level fingerprints, and provided a new discussion on the impact of elastic and viscoelastic rebound on geodetic GPS signals in Greenland. Our results can be summarized as follows.
The long-wavelength (ℓmax= 128) ER of Greenland in response to present-day mass loss shows a clearly bimodal pattern, according to the spatial pattern of the ICESat-derived mass balance model employed in this study (Sørensen et al. 2011). While the rates of vertical uplift at the GrIS margins amount to several millimetres per year with peak values of ∼10 mm yr-’1, in the bulk of Greenland they range between 2 and 4 mm yr-’1. At these wavelengths, we have found no sign of broad scale subsidence in the central portion of the GrIS in response to the mass accumulation evidenced by our preferred M3 mass budget. A direct comparison with results obtained by the GRACE mass balance published by Wouters et al. (2008), performed in the same range of wavelengths (ℓmax= 60) has revealed that the ICESat uplift rates exceed, at the GrIS spatial scale, those of GRACE primarily as a consequence of the largest average mass balance in the time period 2003–2008. At smaller spatial scale, ICESat and GRACE-derived vertical movements mainly differ where mass loss is, according to the M3 ICESat mass balance, spatially more concentrated (this occurs especially in the SG and JI areas, where the misfit between the two uplift maps varies between 2 and 3 mm yr-’1). High-resolution uplift maps based on the mass balance M3 reveal fine details of the elastic uplift pattern of Greenland and of the geoidal variations, including regions that are subject to subsidence in response to local ice re-advance. While in the case of the northeastern Storstrømmen glacier subsidence is expected since the glacier is recovering from a previous surge (Reeh et al. 1994), for Flade Isblink the ICESat record elaborated by Sørensen et al. (2011) possibly constitutes the first evidence of accretion over a multiannual time scale.
The sea level equation provides the natural means of evaluating the contribution of ER to the regional and global sea level variations caused by the melting of the GrIS. We have found that the regional sea level fall along the coasts of Greenland is strongly anticorrelated with vertical movements. In this regional context, where we have confirmed that self-gravitation of the oceans plays a minor role, the approximate equation dS/dt≃-’dU/dt can be applied (Spada et al. 2006). The eustatic (i.e. globally averaged) rate of sea level change associated with the ICESat melting scenario M3 amounts to ≈+0.67 mm yr-’1, in good agreement with the value suggested by the GRACE mass balance recently obtained by Velicogna & Wahr (2006) and Wouters et al. (2008) and also with previous estimates of secular sea level variations based on classical (tide-gauge) observations (Mitrovica et al. 2001). The ICESat-based global sea level signature (fingerprint) shows the typical, non-uniform pattern of near-field sea level fall and far-field sea level rise (well above the eustatic value) expected as the effect of the elastic deformation of the Earth and of self-gravitation (Mitrovica et al. 2001). Except in the immediate near field of the GrIS, the spatial heterogeneity of the ICESat mass balance M3 does not produce significant deviations from the fingerprint associated with an eustatically equivalent (but spatially uniform) mass balance. We have found that the current mass loss in Greenland has a significant impact on key areas as the Mediterranean, where it can produce a sea level rise of ∼0.5 mm yr-’1, comparable with the long-term GIA contribution in the bulk of the basin (Tsimplis et al. 2011).
The GIA component of present-day uplift in Greenland has been evaluated using a suite of different scenarios for the global history of melting since the Last Glacial Maximum, based on different geophysical constraints (Fleming & Lambeck 2004; Peltier 2004). Our findings support previous observations by Velicogna & Wahr (2006) and Khan et al. (2007, 2008), who have emphasized the importance of the GIA component. However, we have uncovered significant breadth in the GIA predictions which is ultimately associated with uncertainties in the ice melting scenarios and mantle rheological profiles. Using combined (ER+GIA) models, in the final part of this work we have addressed the problem of the interpretation of trends of vertical displacement recorded at the five GPS sites with sufficiently long records previously considered in the literature (particularly, by Khan et al. 2007, 2008). We have explicitly shown that ER predictions at these GPS sites are significantly affected by mass loss in near-field drainage basins resolved in the ICESat mass balance and previously excluded in ER modelling (in this respect, a particularly interesting case is that of THUL, where previous ER computations by Khan et al. 2008, have clearly underestimated the uplift rate). According to our computations, only data from two out of the five GPS sites considered (namely THUL and KELY) can be currently employed to constrain competing GIA models; the others (SCOR, KULU and QAQ1) are mainly responding to ER. In the absence of significant tectonic deformations at the GPS sites, we would expect that ER predictions combined with GIA could acceptably explain the observations. According to our modelling efforts, at the sites of THUL, KELY, and SCOR and at the tide gauge of NUUK the GIA (ICE–5G) component of the total rate of displacement effectively helps to reconcile ER (ICESat-based) predictions with the available GPS observations for the time period 2003–2008.
We are grateful to Shfaqat Abbas Khan for kindly providing updated GPS solutions and for valuable discussion. We thanks Erik Ivins and an anonymous reviewer for very constructive suggestions which have helped to improve the manuscript. We are grateful to Valentina Barletta, Spina Cianetti and Matt King for helpful comments, and to Jerry Mitrovica and Thomas James for discussions about the sea level equation. We thank Kurt Lambeck for providing the Australian National University data for the history of continental ice sheets since the Last Glacial Maximum, Bert Wouters for making available the GRACE mass balance and the APLO service for providing the elastic load-deformation coefficients. All the figures have been drawn using the GMT public domain software (Wessel & Smith 1998). The computations involving spherical harmonics have been performed by means of program SHTOOLS http://www.ipgp.fr/wieczor/SHTOOLS/SHTOOLS.html). The numerical codes used for the GIA computations are available from the web (see http://www.fis.uniurb.it/spada/SOFT.html) or from GS. This work was partly supported by COST Action ES0701 ‘Improved Constraints on Models of Glacial Isostatic Adjustment’. The numerical computations have been partly performed thanks to a CINECA award under the ISCRA (Italian SuperComputing Resource Allocation) initiative, for the availability of high performance computing resources and support through the Class C Projects contract n. HP10CS9J50. We acknowledges the ice2sea project, funded by the European Commission′s 7th Framework Programme through grant number 226375 (ice2sea manuscript no. 036).