Present-day mass changes for the Greenland ice sheet and their interaction with bedrock adjustment

Since the launch in 2002 of the Gravity Recovery and Climate Experiment (GRACE) satellites, several estimates of the mass balance of the Greenland Ice Sheet (GrIS) have been produced. To obtain ice mass changes estimates, data need to be corrected for the effect of deformation changes of the Earth's crust. This is usually done by independently modeling the Glaciological Isostatic Adjustment (GIA) trend and then by removing it from the data. Recently, Wu et al. (2010) proposed a new method to simultaneously estimate GIA and the present-day ice mass change, reporting an ice mass loss of around half of the previously published estimates and a general bedrock subsidence concentrated in the central parts of Greenland. This subsidence appears to be counterintuitive since the ice sheet is loosing mass at present. It was suggested by the authors that this could be a new evidence for additional net past ice accumulation. In this study, a 3-D ice-sheet model with a surface mass balance forcing based on a mass balance gradient approach has been used to: (a) analyze the bedrock response to changes in the ice load in order to evaluate whether bedrock subsidence and ice thinning can exist simultaneously; (b) study the magnitude and the pattern of the bedrock movement; and (c) evaluate if present-day bedrock subsidence could be the result of a net past mass accumulation. Under a sine forcing of the annual temperature, that mimics the temperature variations in the Holocene, mass changes yield a delay of the bedrock response of 200 years. Thinning of the ice as well as bedrock subsidence coexist during this period with an order of magnitude equal to the observations by Wu et al. (2010). Although, the resulting pattern of bedrock changes differs considerable: instead of the general bedrock subsidence reported before, we found areas of bedrock uplift as well as areas of bedrock subsidence. A simulation since the last glacial maximum (with the temperature represented as a linear increase from −10 K to present-day) yields a time lag of 1990 years for the bedrock response relative to the temperature forcing and an average uplift of 0.3 mm yr−1 for present-day. The spatial pattern of bedrock-change shows subsidence in the south and northwest as well as uplift in the center and northeast. We obtained these results assuming that the solid earth is a flat elastic lithosphere resting over a viscous relaxed asthenosphere (ELRA model). Using a more sophisticated Self Gravitational Viscoelastic (SGVE) model, we obtain qualitatively similar results: a 2200 years lag and an average uplift for present-day of 0.2 mm yr−1. The spatial pattern of bedrock movement is similar as well. Finally, results are shown for a temperature reconstruction based on ice core data confirming the deglaciation experiment. According to this study, a bedrock subsidence with a maximum in the central parts of Greenland cannot be that recent explained by a net past ice accumulation. This undermines results suggesting that recent loss is only half of the regular ice mass loss changes.

Present-day mass changes for the Greenland ice sheet and their interaction with bedrock adjustment M. Olaizola1,* , R. S. W. van de Wal 1 , M. M. Helsen 1 , and B. de Boer 1

Introduction
Since the Fourth Assessment Report by the International Panel of Climate Change, which included a mass balance (MB) estimate for the Greenland ice sheet (GrIS) of −50 to −100 Gt yr −1 from 1993 to 2003 (Lemke et al., 2007), several studies have been published showing an increased loss in MB.Estimates mostly range between −160 to −300 Gt yr −1 over the period 2002 to 2010 (Velicogna and Wahr, 2006;Wouters et al., 2008;Rignot et al., 2008;Velicogna, 2009;Tedesco et al., 2010).These values are in agreement with the study presented by van den Broeke et al. (2009), where mass balance estimates derived from Gravity Recovery and Climate Experiment (GRACE) data are compared with the mass budget method.
Recently, a lower MB estimate (−81 to −127 Gt yr −1 over the period 2002 to 2008) inferred from GRACE data has been reported by Wu et al. (2010) (referenced as W10 hereafter).The authors proposed a new method to estimate simultaneously the Glaciological Isostatic Adjustment (GIA) trend and the present mass loss from GRACE, while in previous studies the GIA trend is removed a priori from the data (e.g., Velicogna and Wahr, 2006;Velicogna, 2009).Using this new method, W10 obtained that thinning of Introduction

Conclusions References
Tables Figures

Back Close
Full ice at the margin of the GrIS is stronger and more extended than the thickening at the center, in agreement with previous results (Thomas et al., 2006).The authors also retrieved bedrock subsidence almost over the entire island, most pronounced in the center and with a small area of bedrock uplift in the northeast.Averaged over Greenland, the bedrock displacement is estimated at −0.56 ± 0.17 mm yr −1 .In this study we tested the hypothesis that the recent negative trend in the integrated SMB over the GrIS intuitively leads to ice thinning and hence an average uplift of the bedrock response, which is in disagreement with the results by W10.Therefore, we analyze whether bedrock subsidence and ice thinning can exist simultaneously for the presentday (PD) Greenland Ice Sheet, and whether the sign, magnitude and pattern of the vertical bedrock movement as reported by W10 can be simulated by realistic reconstructions of the ice-sheet evolution.This is done with a coupled ice sheet-bedrock model driven by variations in mass balance.
We start with a description of the three different components of the 3-D ice sheetbedrock model: the ice dynamics, the surface mass balance and the solid earth.In the following section we present the results of different simulations where we studied the relation between changes in ice thickness and the bedrock response.The first experiment schematically mimics climate fluctuations during the Holocene following a sine function, and describes a quasi steady-state behavior.To incorporate the influence of the last glacial era, we carried out a second experiment adding a linear increase in temperature from glacial to PD conditions.Those results are performed for two different bedrock models.Finally, we implemented a temperature forcing based on ice core records for a more realistic experiment.Introduction

Conclusions References
Tables Figures

Back Close
Full and atmosphere are treated with the "SMB Gradient Method" (Helsen et al., 2011), that allows the use of SMB fields from climate models and accounts for the height-mass balance feedback.The SMB from the regional climate model RACMO 2/GR (Ettema et al., 2009) is used as the reference SMB field, and perturbations are applied based on local changes in surface elevation (H s ) that results from the ice-sheet model.To enable an additional perturbation due to external temperature forcing, variations in the annual surface temperature ∆T are divided by γ, the temperature lapse rate that results from a linear regression between elevation and temperature (γ = −7.3743K km −1 ).In this way, ∆T is translated into a change in surface elevation (∆H = ∆T/γ).Taking this into account, the SMB is calculated by: where a and b are local coefficients.In this method there is a distinction between the accumulation regime and the ablation regime.Therefore, two regression equations are needed (two b and two a coefficients), which are obtained by minimizing the perpendicular distances between data points and regression lines (Helsen et al., 2011).
In Fig. 1, the two best-fitting curves for a point located at the west coast are shown.
To get the coefficients a and b for a certain grid point in the GrIS, it is necessary to consider the SMB, determined by RACMO, of all the neighboring points in a radius of 150 km (Helsen et al., 2011).To avoid the positive SMB to increase or decrease without restrictions, maximum and minimum values are calculated (Helsen et al., 2011).The For the second component (the ice-sheet), ice dynamics and thermodynamics are treated solving the equation for mass conservation (van der Veen, 1999): where H i is the ice thickness and Ū the vertically averaged horizontal velocity.The temperature field is solved by: where k and c p are the thermal diffusivity and specific heat capacity of the ice.The dissipation of energy due to ice deformation is represented by Φ.If the basal temperature reaches the pressure melting point, ice sliding over the bedrock is considered with a Weertman-type sliding law (Weertman, 1964).
The interaction with the the solid earth accounts for the bedrock response to the loading changes of the overlying ice.Changes in bedrock height modify the ice sheet surface elevation and therefore the surface mass balance and ice thickness.We use two different formulations for the bedrock response.In the first approach, the Earth is assumed to be a flat elastic lithosphere (EL) resting over a viscous relaxed asthenosphere (RA).This is the so called ELRA model (Le Meur and Huybrechts, 1996).Accordingly to the ELRA model, the bedrock responds with a downward deflection w 0 to the pressure exerted by a point load q.The steady state displacement w 0 is given by the following equation for a normalized distance x = r/L r from q (Le Meur and Huybrechts, 1996): where D is the flexural rigidity (1 × 10 25 Nm, a typical value for land) that allows the depression to extend beyond the point q where the load is located, χ is the zero order Introduction

Conclusions References
Tables Figures

Back Close
Full .The elasticity of the lithosphere is assumed to be linear, hence, the total deflection of the bedrock at some point can be calculated as the sum of the contribution of all neighboring points.While the lithosphere accounts for the shape of the deformation, the asthenosphere controls the time response.The rate of the vertical bedrock movement, ∂w 0 ∂t is proportional to the deviation of the profile from the equilibrium state, w − w 0 , and inversely proportional to the relaxation time τ (Le Meur and Huybrechts, 1996): The value of τ is fixed to 3000 years.
The second approach is a more complete physical approach, the Self Gravitational Viscoelastic (SGVE) model, consisting of an elastic lithosphere, two viscoelastic mantle layers and an inviscid core (Le Meur and Huybrechts, 1996).In this model, the isostatic response to a load L s (x,y,t) is given by Y (x i ,y j ,t) where ∆x, ∆y are the spatial resolution, ∆θ is the angular distance between grid points and T m is the preceding 20 000 years of ice-loading history.We use the SGVE model in one experiment to compare and validate results from the default ELRA runs.
To initialize the model, bedrock topography and ice thickness were taken from Bamber et al.In order to analyze if a current bedrock subsidence is a possible response to the changes in ice thickness, we carried out a series of experiments, varying the forcing, the simulation duration, and the way the solid earth is modelled.

Last millennium
In the first experiment, we mimicked temperature variations in the last millennium by a sinusoidal function that oscillates around zero with an amplitude of 1 K and a period of 1000 years.This approximation is in reasonable agreement with the reconstruction from Kobashi et al. (2009).The evolution of the mean time derivatives of ice thickness ( H i ) and bedrock elevation ( H b ) are shown in Fig. 2a (spatial average of the variables).
We plot the time derivatives to illustrate more precisely the phase lag between ice thickening (thinning) and bedrock subsidence (uplift).The temperature forcing (indicated without units by the black line as a reference) leads to ice thickness changes (green), which in turn lead to a bedrock response (blue).We focused the analysis on the last period of the run, once the system reaches a new quasi steady-state after 57 kyr.The total length of the simulation is 60 kyr.
In periods when the temperature is above zero, H i has negative values, corresponding to ice thinning and the bedrock responds with uplift (positive values of H b ) after some time.For this experiment there is a time lag of 200 yrs, represented by the grey area in Fig. 2a.During this time period, ice thinning and bedrock subsidence exist simultaneously.The time lag and its possible variations, as a result of changes in the forcing or in the way the solid earth is modelled, are important to consider as it might be an explanation for the results by W10.
For this experiment we obtained a maximum average bedrock subsidence of H b = 0.4 mm yr −1 .This is the spatial average over Greenland and higher values can be found near the margins as shown in Fig. 2b, c, d and e, where we present the spatial patterns of H b .Figures 2b and c  to the time slices where H i reaches its maximum and minimum values respectively (green points in Fig. 2a).In Fig. 2d and e (t=59 100 yrs and t=59 400 yrs), H b is at its minimum and at its maximum respectively (blue points in Fig. 2a).
The pattern of H b at t=57 900 yrs (Fig. 2b) is characterized by bedrock subsidence along the western margin.In this ablation area, the SMB increases as a consequence of lower temperatures.As a result, ice thickens and the bedrock subsides (blue colors).In the southeast, an accumulation area, the opposite occurs: a decrease in precipitation due to lower temperatures is translated into thinning of the ice and therefore into bedrock uplift (yellow to red colors).When the temperature increases to positive values at t=58 250 yrs, the bedrock subsides in the accumulation area located in the southeast.This is due to the enhanced precipitation that results in ice thickening.Over the rest of the island, the temperature increase is followed by an increase in ablation, which results in less ice volume (minimum value of H i ) and bedrock uplift.At t=59 100 yrs the pattern is similar to the one at t=57 900 yrs.Bedrock subsidence is observed along the western margin because the bedrock is still reacting to past changes when the temperature was below zero.This results in an ablation reduction, thus ice thickening.In the southeast, precipitation decreases, ice thins, and the bedrock reacts with uplift.At t=59 400 yrs, H b reaches its maximum and the opposite occurs.The pattern is similar to the one at t=58 250 yrs, characterized by bedrock uplift in the southwest and bedrock subsidence in the southeast.
In all cases, the pattern of the signal is most pronounced in the southern half of Greenland.Additionally, the signal alternates periodically between bedrock subsidence and uplift being stronger along the margin.This is due to the spatial pattern of surface mass balance with an ablation area located in the south west and an accumulation region in the southeast.As such results do not depend much on the details of the mass balance formulation.Introduction

Conclusions References
Tables Figures

Back Close
Full

Holocene and last deglaciation
Since bedrock response shows a lag to changes in ice thickness, it is also interesting to study the effect of the last glacial era on the magnitude of the lag.In order to do this, we carried out a steady-state simulation with PD temperatures lowered by 10 K for 100 kyr.
After this, we forced the model with a linear increasing temperature to PD in 10 kyr.
Finally, the temperature oscillates around PD as a sine function with an amplitude of 1 K and a period of 1000 years for another 10 kyr, similar to the previous experiment.
In Fig. 3a the forcing is illustrated as a black line, the green line shows changes in ice thickness and the bedrock response is shown in blue.
During the first 3000 years (between −20 000 years and −17 000 years), the temperature increase results in more ablation and, hence, ice thinning (negative values of H i ).This is a consequence of more ablation.Hereafter, warmer conditions (between −4 K and PD temperature) also enhance precipitation, which becomes the dominant effect on this period.Therefore, ice starts to thicken (positive values of H i ) around the year −17 000.After 9000 years (left side of the grey rectangle in Fig. 3a), H i changes sign to negative values, meaning that ice thins due to the increasing temperature.The bedrock reacts to this ice load reduction with the expected uplift, although with a delay of 1990 yrs, indicated by the grey rectangle in Fig. 3a.Hence, there is a period where ice thinning coincides with bedrock subsidence.The spatial pattern of H b at the moment of maximum average subsidence, at t=−11 000 yrs (Fig. 3b) is characterized by subsidence almost on the entire island, with higher values in the center.Some areas with uplift are present as well, mainly along the western ablation zone.The average value at this time is H b = −0.9mm yr −1 .The strong signal of bedrock subsidence in the center of Greenland diminishes with time as can be seen in the spatial patterns of bedrock elevation presented for the years −8800 Since different physical formulations of the solid earth might have a strong impact on the values of H b , we used the more sophisticated Self Gravitational Viscoelastic (SGVE) model in addition to the ELRA model, to compare and validate the previous results.
Figure 4a shows the values of H i (left plot) and H b (right plot) for the two models.In the initial period the results with the SGVE model show some unrealistic values, which are due to a lack of spin-up time of the SGVE model: this model has a memory for the previous 20 kyr but at the beginning of the simulation, the starting file is the result of a previous steady-state run with a temperature 10 K lower than PD temperature, and the past memory is lost.Therefore the model needs a spin up period, after which the results with the SGVE model are similar to the ones with the ELRA model.The general behavior of the H b patterns (Figs. 3 and 4) are similar as well: as a result of the last glaciation there is a strong bedrock subsidence in the center that vanishes with time.For the SGVE model, the maximum value of H b = −0.8mm yr −1 occurs at t=−11 000 yrs (similar to the ELRA model), and then starts to diminish, reaching slightly positive values around t=−8800 yrs with H b = 0.05 mm yr −1 .At t=−5000 yrs, the maximum bedrock uplift takes place, H b = 0.3 mm yr −1 , after which the average value decreases again to H b = 0.2 mm yr −1 for PD.The largest differences in the patterns for the two models occur in North Greenland.Bedrock uplift is present in this region for the last three selected time slices in the SGVE model, while with the ELRA model, there is practically not bedrock movement in that area.Moreover, the time lag of the bedrock response increases slightly with respect to the ELRA model to a value of 2.200 yrs.

Holocene and last deglaciation, ice core data
To study the bedrock response to PD changes in ice thickness in a more realistic way than in the previous schematic experiment, we used as forcing, a temperature record from an ice core based on the GRIP δ 18 O converted into a surface temperature deviation following Johnsen et al. (1995).Sea level reconstruction is prescribed from Introduction

Conclusions References
Tables Figures

Back Close
Full deconvoluted marine benthic oxygen isotopes from Bintanja and van de Wal (2008).Results as obtained using the ELRA model are shown in Fig. 5a.Between the years −20 000 and −10 000 temperatures remain below zero, which allows ice thickening, especially when the temperature strongly increases (as for t=−14 500 yrs), because a warmer climate is characterized by higher precipitation rates.In fact, for very cold periods (from −20 000 to −15 000 years), precipitation is reduced, which yields ice thinning and bedrock uplift.During the Holocene (from t=−10 000 yrs to PD), temperature oscillates around zero.This results in ice thinning from t=−10 200 yrs (indicated by the left side of the grey rectangle in Fig. 5a).Negative values of H i persist during the Holocene with a few short periods of ice thickening.The ice load reduction causes bedrock uplift (t=−8400 yrs, right side of the grey rectangle) that occurs with a lag of 2800 yrs.
In Fig. 5b, c and d

Conclusions
To analyze the interaction between variations on the ice load and on the bedrock response in the Holocene, we carried out the experiment presented in Sect.3.1.We found a time lag of the bedrock response of 200 years, implying that the bedrock might be reacting to a past accumulation change that coincides with the Little Ice Age, between the years of 1400 and 1900.The time lag of the bedrock response allows for a situation of simultaneous ice thinning (due to an increase in temperature or negative MB) and bedrock subsidence as suggested by W10.

TCD Introduction Conclusions References
Tables Figures

Back Close
Full We obtained similar values for this lag varying the amplitude (up to 5 K), as well as varying the bedrock relaxation time τ to 5000 yrs, which was fixed to 3000 yrs in the experiments presented in this paper.In the ELRA model, the long term response of the bedrock is controlled by the asthenosphere, being inversely proportional to the relaxation time.The numerical values of τ are based on assumptions of the viscosity of the Earth interior.Moreover, about half of the error in the models used to calculated GIA trends is due to the lack of information of the earth viscosity profile (Velicogna and Wahr, 2006).Variations in τ introduce a maximum difference in the values of the lag of 20 % (in case of a period in the forcing of 10 000 years).Therefore, an error of 20 % can be assigned to the magnitude of the time lag.
The maximum average bedrock subsidence found in this first experiment is H b = −0.4mm yr −1 , increasing up to −1 mm yr −1 near the margins.Although we can obtain a result for which the order of magnitude of the subsidence is in accordance with the one reported by W10, the spatial pattern of bedrock changes differs considerably.We observed the highest values in the lower half part of the island and not in the center.
Moreover, oscillations between subsidence and uplift in the southwest and southeast occur instead of a general subsidence, because in those areas the largest ice changes take place.In fact, any realistic mass balance forcing will lead to a stronger response along the margin than in the center as ice thickness changes are larger near the margin.
Nevertheless, it should be noted that this first experiment is a schematic way to approximate the present-day conditions of the GrIS.A limitation of this experiment is the steady-state initial condition.Therefore a second experiment was carried out with the last (de)glaciation preceding the Holocene temperature variations (Sect.3.2).Also in this case, a period of time exists where bedrock subsidence and ice thinning occur simultaneously.Furthermore, the spatial pattern of bedrock elevation shows a strong subsidence in the center just after the deglaciation started, even though in combination with uplift in the west, which is an ablation area.The intensity and the extension of the central subsidence decrease over time until they almost vanish.For PD conditions, Introduction

Conclusions References
Tables Figures

Back Close
Full In Table 1, we show a summary of the results of the ice sheet-bedrock experiments presented in this paper.In all the simulations we applied the novel SMB gradient parameterization (Helsen et al., 2011), where the time evolution of the SMB is based on changes in surface elevation rather than via a constant lapse rate as is the case for the PDD method, often used in ice-sheet models.The SMB model formulation has an influence on the results.This is clear in the first experiment (Sect. 3.1) where an alternation between bedrock subsidence and uplift in the southwest and southeast of Greenland was found.This is due to the SMB field which is characterized by strong ablation in the southwest and an accumulation area in the southeast.
A limitation of the applied model is the lack of detail with respect to outlet glaciers, particularly in the southeast.Those outlet glaciers are partly in different regions, so the thickness change pattern maybe somewhat different from what we find, but the largest changes will still take place in the marginal zones and yield qualitatively similar results for the bedrock response pattern.
In all the experiments we found a time lag of the bedrock response not higher than 3000 years.This implies that for the PD conditions, after 10 kyr of deglaciation, the bedrock is adjusted to the ice load reduction and an average bedrock uplift is present in Greenland.Therefore, we conclude that for present day, a bedrock subsidence with a maximum in the central parts of Greenland as reported by W10 cannot be explained by past mass changes in the surface mass balance as the authors suggested, not by the deglaciation, but also not by changes since the Little Ice Age.This undermines the result of a mass loss of half of the values reported in earlier studies.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | transition from one regime to the other is defined by a critical surface elevation H c , obtained when the SMB for the accumulation and the ablation regimes are equal.As the simulation evolves over time, the initial altitude of a grid point is modified by changes in bedrock topography, ice thickness, or by changes in surface temperature (translated into changes in surface elevation via Eq.1).When surface elevation crosses H c , a switch of regimes occurs and therefore a switch of the values a and b used to calculated the SMB.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Kelvin function, r is the original distance from the load point and L r the radius of relative stiffness: L r = D/ρ a g 1 4 (2001)  with a spatial resolution of 20 km × 20 km.
(t=57 900 yrs and t=58 250 yrs) correspond Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and −5000, with values of H b = −0.1 mm yr −1 and H b = 0.4 mm yr −1 respectively.For PD conditions (t = 0 in this experiment), the bedrock subsidence at the center almost vanishes and the pattern is characterized by an average uplift, H b = 0.3 mm yr −1 .Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | , we show the patterns of bedrock elevation changes for the selected red points in Fig. 5a.The maximum bedrock subsidence is at t=−9980 yrs, reaching a value of H b = −1.3mm yr −1 .As in the previous more schematic experiments, the strong bedrock depression in the center of Greenland fades with time.At t=−8000 yrs the bedrock subsidence is reduced to H b =-0.2 mm yr −1 .Later on, at t=−2000 yrs, the maximum bedrock uplift occurs with a value H b = 0.8 mm yr −1 , which decreases to H b = 0.5 mm yr −1 for PD conditions.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

'Hb´Fig. 3 .
Fig. 2. (a) Time evolution of the spatial average H b (blue line) and H i (green line) for a temperature forcing (black line) simulating in a schematic way the temperature variations at the Holocene.We present the last cycles were a new steady-state is reached.As temperatures increase, ice thins (negative values of H i ) and the bedrock reacts with uplift with a lag of 200 yrs indicated by the grey area.The spatial patterns of bedrock elevation are presented for a moment where: (b) H i reaches its maximum value; (c) H i is at its a minimum, (d) H b is at its minimum, (e) H b is at its maximum.
Fig. 5. (a) Time evolution of the spatial average H b (blue line) and H i (green line) for a temperature forcing from ice core data (black line).The grey rectangle shows the time lag of the bedrock response of 2,800 yrs.For the selected red points, the spatial patterns of bedrock elevation are presented for (b) the maximum subsidence; (c) a moment when ice thinning and bedrock subsidence exist simultaneously; (d) a moment close to the maximum uplift; (e) present-day conditions, characterized by an average bedrock uplift.

Table 1 .
Peltier (2004) and lag values as well as H b and the PD-bedrock pattern obtained in this study and the ones reported byPeltier (2004) and Wu et al. (2010).