Interactive comment on “ Changing basal conditions during the speed-up of Jakobshavn Isbræ , Greenland ” by M . Habermann

Abstract. Ice-sheet outlet glaciers can undergo dynamic changes such as the rapid speed-up of Jakobshavn Isbrae following the disintegration of its floating ice tongue. These changes are associated with stress changes on the boundary of the ice mass. We invert for basal conditions from surface velocity data throughout a well-observed period of rapid change and evaluate parameterizations currently used in ice-sheet models. A Tikhonov inverse method with a shallow-shelf approximation forward model is used for diagnostic inversions for the years 1985, 2000, 2005, 2006 and 2008. Our ice-softness, model norm, and regularization parameter choices are justified using the data-model misfit metric and the L curve method. The sensitivity of the inversion results to these parameter choices is explored. We find a lowering of effective basal yield stress in the first 7 km upstream from the 2008 grounding line and no significant changes higher upstream. The temporal evolution in the fast flow area is in broad agreement with a Mohr–Coulomb parameterization of basal shear stress, but with a till friction angle much lower than has been measured for till samples. The lowering of effective basal yield stress is significant within the uncertainties of the inversion, but it cannot be ruled out that there are other significant contributors to the acceleration of the glacier.


Introduction
Ice-sheet outlet glaciers can evolve much more dynamically than formerly thought (Truffer and Fahnestock, 2007).Modeling and understanding the processes involved in these rapid changes is challenging.Despite the abundant surface data available from satellites, conditions within the ice and at the base of the ice are still difficult to observe, but these are crucial components of successful prognostic ice-sheet models.
Jakobshavn Isbrae is one of the most active outlet glaciers in Greenland and has a century-long record of observations (Weidick et al., 1990).This outlet glacier drains about 5.5 % of the ice-sheet area (Rignot and Kanagaratnam, 2006) and has undergone a rapid evolution in the last two decades.During the 1990s Jakobshavn Isbrae had a relatively stationary terminus position (Sohn et al., 1998), but starting in 1997, increased thinning of the floating ice tongue was observed (Thomas et al., 2003), followed by the retreat and complete disintegration of the 15 km-long ice tongue in 2003 (Podlech and Weidick, 2004).Coinciding with the retreat of the ice front, the ice underwent a significant speed-up, almost doubling its speed by 2003 (Joughin et al., 2004).After the disintegration of the ice tongue, the ice front retreat and the accelerations in speed have decreased but are still ongoing today (Joughin et al., 2012).
Three main processes have been identified that can contribute to the changes in outlet glaciers generally and at Jakobshavn Isbrae specifically (Joughin et al., 2012).The first process is a speed-up of the ice to compensate for a loss of buttressing (downstream contact with the bed and/or fjord walls) during the retreat of the ice front.The relationship between front position and speed has been well observed on longer timescales and on seasonal timescales (Joughin et al., 2008b;Amundson et al., 2010).The second process is a loss of overburden pressure through thinning of the ice, while the basal water pressure is assumed to be fixed through its connection to the ocean.This leads to a decrease in effective pressure and a decrease in basal shear stress, which in turn leads to an increase in sliding speed (Meier and Post, 1987;Pfeffer, 2007).The third process is a steepening of slopes induced by the strong thinning on the main trunk, causing the speed-up to diffuse inland (Joughin et al., 2008b;Payne et al., 2004).Other possible processes include weakening of the ice in the lateral shear margins and increase in basal water pressure through changes in the hydrological system (Van Der Veen et al., 2011).The observational evidence strongly favors an acceleration mechanism that is ocean and terminus driven (Motyka et al., 2011;Joughin et al., 2012).
The well-observed changes of Jakobshavn Isbrae make it possible to investigate temporal changes in effective basal yield stress by inverting surface velocities for different years.(Joughin et al., 2012) performed one inversion for the 1990s velocities and one for the 2009 velocities.Here we expand on this by inverting all available velocity fields and by conducting an extensive parameter study to discuss the robustness of the inversion results.
To take advantage of the wealth of surface data, we use inverse methods to reconstruct conditions at the icebed boundary.Inverse methods were first introduced to the field of glaciology by (MacAyeal, 1992), and have since been used, improved and extended in multiple studies (e.g., Truffer, 2004;Maxwell et al., 2008;Raymond and Gudmundsson, 2009).Much like other recent studies (Morlighem et al., 2010;Konovalov, 2012;Petra et al., 2012), we use a Tikhonov regularization to stabilize the solution, and we focus on justifying the choices that accompany this method.
In this study we investigate different parameter choices for the effective basal yield stress inversion of Jakobshavn Isbrae, where decisions are mostly based on the data-model misfit metric.The chosen parameters are then used to invert, for effective basal yield stress, the surface velocity data sets of the years 1985, 2000, 2005, 2006 and 2008.We discuss the robustness of these results and the agreement with commonly used parameterizations of effective basal yield stress.

Model
To investigate spatial changes and characteristics of basal shear stress, we use the shallow-shelf approximation (SSA) (Morland, 1987) as the forward model in a Tikhonov inversion.

Forward model
The Parallel Ice Sheet Model (PISM) is a 3-D thermomechanically coupled hybrid ice-sheet model that solves a combination of the shallow-ice and shallow-shelf approximations (Bueler and Brown, 2009, http://www.pism-docs.org).In this study only the SSA is used and the vertically averaged ice softness does not vary horizontally.Details about the SSA can be found in Schoof and Hindmarsh (2010) and the implementation in PISM is described in Bueler and Brown (2009).
We follow Joughin et al. (2012) and use the SSA as a forward model.Despite being depth-averaged the model does consider membrane stresses, vertical shear on the other hand is not considered.Ignoring vertical shear can be justified by the weak temperate basal ice layer that is present at Jakobshavn Isbrae, which concentrates vertical motion near the bottom, and by the weak bed compared to the driving stresses, which leads to motion that is dominated by basal ice motion, at least in the lower regions of the glacier (Lüthi et al., 2002).However, it is important to keep in mind that the results derived in this paper are effective basal yield stress fields that are consistent with the SSA and surface observations, and might not reflect actual physical till properties.
The input fields needed for the forward SSA are ice thickness H , surface elevation z s , ice softness A, and a basal shear stress τ b .The model output is the surface velocity u.PISM treats the SSA as if it applies to the entire grid domain, even in ice-free locations.Each grid point can be either icy or ice free, and either grounded or ocean, for a total of four states.A point is ice free if the ice thickness H falls below a small threshold (set to 0.01 m).The distinction between ground and ocean is made by computing what the surface elevation would be at that location for grounded ice and for floating ice; the maximum elevation determines the state.In regions where H is zero, the product of effective viscosity and thickness is regularized with a constant (set to 1 × 10 13 Pa s m), for details see http://www.pism-docs.org.The value of τ b is adjusted based on the ice/ice-free grounded/ocean status of a grid point.For floating locations, the value is set to 0, and for ice-free ground it is a large constant.Consequently, τ b depends on the effective yield stress τ c only for grounded ice.Dirichlet boundary conditions (i.e., locations where υ is known) are prescribed at the outermost 5 km of the rectangular domain.To approximate the ocean front boundary condition for the shallow-shelf approximation, we extend a thin notional shelf into all ice-free areas (see (Bueler and Brown, 2009, Sect. 2.6) for details).No additional boundary conditions are applied to the terminus of the glacier, instead the ice thickness simply decreases to zero from one grid point to the next.In this way the glacier outline is determined by the ice thickness given in the DEM for each year.The change in buttressing forces is implicitly taken into account by adjusting the ice geometry.In this paper we use these fixed approximations of the buttressing forces and we invert for different distributions of effective basal yield stress.We chose a grid resolution of 500 m × 500 m.A finer resolution is not warranted by the data and tests with coarser grids show convergence.A finer grid might be desirable in the area of the deep trough, where basal topography changes rapidly.
The basal shear stress τ b is parametrized through a power law: where u is the basal sliding velocity, and the threshold velocity u threshold is set to 100 m a −1 .The purely plastic case is achieved by setting q = 0, whereas q = 1 leads to the common treatment of basal till as a linearly viscous material: τ b,x = γ u and τ b,y = γ v, where γ ≥ 0 is a scalar function of position, called the basal stickiness.When setting q = 1 the basal stickiness, γ , and the effective basal yield stress, τ c , are related through γ = τ c u threshold .Here, instead of setting q = 1 and solving for γ we solve for τ c , which has units of stress and is the basal yield stress if q = 0. We approximate the perfectly plastic case by setting q = 0.25 for this study and call τ c the effective basal yield stress.This is the stress that occurs at the base when it is sliding at u threshold .If q is close to zero, and if there is reasonable sliding at the base, we can expect that the stress at the base will be close to τ c .As such it plays the role of a yield stress.Test inversions with q = 0.1 and q = 0.001 for the 1985 and 2006 data sets result in different τ c values, but the pattern and amplitude of changes in τ c remain and the main conclusions of this paper are unchanged.The positivity of τ c is enforced by solving for ζ in τ c = τ c,scale exp(ζ ) where τ c,scale is a scale parameter to keep ζ of order 1 for typical values of τ c .
The chosen values for q and u threshold used here were found to provide the best representation of observed ice motion (Bueler, personal communication, 2012).As mentioned before, the results derived in this paper are effective basal yield stress fields that are consistent with our model choices and surface observations, and might not reflect actual physical till properties.The main conclusions of this paper, namely a weakening of the till near the terminus, remain valid for different choices of q and u threshold .
We assume the instantaneous (diagnostic) surface velocities represent instantaneous deformation rates and effective basal yield stress at depth.In other words, no timedependent (prognostic) runs are performed and instead the forward model calculates a velocity field from effective basal yield stress τ c , and the inversion is an attempt to recover τ c from measured surface velocities at a given time.

Inferring effective basal yield stress
Solving for the effective basal yield stress distribution is an ill-posed inverse problem, one consequence being the multitude of possible solutions.Often these ill-posed problems can be stabilized by imposing additional constraints that bias the solution.This is referred to as regularization (Aster et al., 2005).We apply the widely used Tikhonov regularization, which defines a cost functional, I (τ c , α), with an added reg-ularization term: where M is the data-model misfit, N is the model norm (regularization term) and α is the regularization parameter.Note that, depending on the application, α is sometimes attached to the model norm instead of the data-model misfit.This only changes the value of α, but not any of the results.We discretize the functional I (τ c , α) by representing τ c via a finite-element approximation, and by computing a finite element solution for u(τ c ). Doing so determines a discretized functional I disc : R n → R, where n is the number of grid points where τ c is defined.The gradient of this discretized functional (with respect to the standard inner product on R n ) can be computed exactly; the gradient of the term M 2 is determined by a discrete computation similar to the continuous computation of the L 2 gradient described in the Appendix of (Habermann et al., 2012), whereas the gradient of N 2 (which is quadratic in τ c ) is computed trivially.A minimum of the discrete functional can then be sought by any one of a number of gradient-based minimization algorithms.We use a limited-memory, variable-metric method from the Toolkit for Advanced Optimization (TAO) (Munson et al., 2012) to seek an exact minimum of the discretized cost function, I (τ c , α).
Assuming that there is a unique minimum (which is true at the very least when α is small), an exactly computed minimum of the discretized functional will be independent of the numerical method used to find it.The area is defined by grounded ice (determined by hydrostatic equilibrium) and the consistent availability of velocity observations over the time periods considered.This is only part of the model domain (see Fig. 1), but all interpretations will be restricted by it.Below we refer to as the 'misfit area'.The model norm in Eq. ( 4) is composed of two parts: the Euclidian L 2 norm and a Sobolov H 1 norm that measures the function's roughness.The factors c L 2 and c H 1 determine the relative weights of these two norms.The variable K defines a typical length scale to rescale the H 1 norm (set to 5 × 10 4 m).The model norm is measured as a difference from a prior estimate τ prior c .A choice of c L 2 = 1 and c H 1 = 0 results in a pure L 2 model norm, which gives preference to solutions with a small departure from the prior estimate.At the other end of the spectrum, setting c L 2 = 0 and c H 1 = 1 results in a pure H 1 model norm, which biases the solution towards smooth differences to the prior estimate.
Achieving a better data-model misfit M carries the cost of a larger model norm.Each choice of the regularization parameter α determines a unique value for the data-model misfit and hence the model norm.To discuss the choice of regularization parameter, α, we introduce the following vocabulary.The observation error is defined as T obs , the system error is defined as T tot = T mod + T obs , where the modeling error, T mod , contains errors from model simplifications and errors in input parameters such as ice geometry.For an ill-posed inverse problem it is not desirable to find an exact minimizer of the data-model misfit, M, because this would lead to overfitting of the data (Habermann et al., 2012).The achieved data-model misfit should not be smaller than the combined error of observations, model simplifications, and parameter choices, T tot .On the other hand, if the data-model misfit is too large, because we are forcing a high degree of smoothness in the effective basal yield stress solution, the highest possible resolution is not achieved and the data are underfit.
There are different ways to choose the regularization parameter α.The "discrepancy principle", which sets the datamodel misfit equal to T tot is useful in situations where all errors in the system are known or where the observation errors can be estimated and the model errors are negligible.For the Tikhonov regularization the discrepancy principle cannot be applied directly.Instead a value for the regularization parameter α is chosen and the resulting data-model misfit value is compared to T tot , if it is known.
A more common situation arises when the errors in the system are not known.It is particularly difficult to quantify model errors that originate from the use of lower order forward models, such as the SSA, and the effect of poorly constrained model parameters, such as the ice softness and bed topography, that are not part of the inversion procedure.In such cases it is possible to use a heuristic "L curve" method (Jay-Allemand et al., 2011;Gillet-Chaulet et al., 2012).It has been proposed for its ease of use, despite some poten-tial shortcomings (discussed in e.g., Vogel, 1987, ch. 7).In the L curve method the data-model misfit is plotted against the model norm (either on a log-log or a linear scale).This curve typically has an L shape and the regularization parameter value corresponding to the "corner" of the curve, which is usually defined as the point of highest curvature, is chosen.The rationale behind this choice of regularization parameter is that past this corner even a small improvement in the datamodel misfit can only be achieved through a large increase in the roughness of the solution.
The actual value of the data-model misfit depends on the misfit area.Therefore, the data-model misfit value can only be used to compare different inversion results if the misfit areas are identical.Here we use the same misfit area for all years, given by the consistent availability of velocity observations and by grounded ice (the 2008 grounding line limits the misfit area in the terminus region).This misfit area is shown in Fig. 1.An appropriate data-model misfit can still lead to overfitting in some subareas and underfitting in others.

Data
A combination of previously published airborne and spaceborne data sets, collected between 1985 and 2008, are used as input to the model.All data sets are given on or interpolated to a 500 m × 500 m grid, which is the grid size chosen for the model.Table 1 gives a summary of surface elevations and velocity fields used for each year.

Surface elevation
We used the 1985 and 2007 digital elevation models (DEM) derived by Motyka et al. (2010).The 1985 DEM is based on aerial photos, whereas the 2007 DEM was derived from SPOT-5 imagery under the SPIRIT (stereoscopic survey of Polar Ice: Reference Images and Topographies) Polar Dali Program (Korona et al., 2009).To extend the model domain, we took lower resolution surface elevations given by Bamber et al. (2001), and substituted the high resolution DEMs in the coverage area.As a result there are sharp transitions from the high resolution DEM to the low resolution DEM.These sharp transitions result in unphysical driving stresses and we smooth the DEM by performing a short (2 week) nonsliding shallow-ice approximation run on a regional scale with PISM.The model domain was chosen beyond the extent of the high resolution DEMs to minimize the impact of boundary effects on the results.Model results are only evaluated within the coverage area of the high resolution DEMs (Fig. 1).
For the years 2000-2008, we used the 2007 DEM together with annual elevation-difference maps from Joughin et al. (2012).
Table 1.Summary of velocity fields and surface elevation data sets used for each year, including details on acquisition dates and source references.The 2007 SPOT DEM that is mentioned was obtained 24 July 2007.The outline of the glacier is given by the ice thickness of the DEM for each inverted year.The misfit area is the same for all years (see Fig. 1).

Year Period covered by vel. field
Reference for vel.Date of surface DEM Reference for DEM

Bed elevation
The bed DEM was developed at the University of Kansas using data collected by their airborne depth-sounding radar (Plummer et al., 2008).It is important to point out that the bed elevation is one of the model input fields with significant uncertainties.Even though the Jakobshavn Isbrae drainage area has been flown repeatedly with a radar depth sounder, the deep trough with its steep margins often does not allow for clear bed returns.
We investigate the influence of bed topography on the inversion results in (Habermann, 2013) and we find that errors in bed topography lead to residuals that are larger than the residuals due to errors in velocity observations.This large expected error is consistent over all inversions performed here and we do not expect a significant influence on the changes in effective basal yield stress.

Ice flow velocity
NASA's Making Earth System Data Records for Use in Research Environments (MEaSUREs) program, provides annual ice-sheet-wide velocity maps for Greenland, derived using Interferometric Synthetic Aperture Radar (InSAR) data from the RADARSAT-1 satellite.The data set contains ice velocity data for the winter of 2000-2001 and 2005-2006, 2006-2007, and 2007-2008 acquired from RADARSAT-1 InSAR data from the Alaska Satellite Facility (ASF), and a 2008-2009 mosaic derived from the Advanced Land Observation Satellite (ALOS) and TerraSAR-X data (Joughin et al., 2010).Here we are using all available velocity data sets except for 2007-2008, which contains data gaps.
For the 1985 inversion we use a velocity data set derived from feature tracking of orthophotos used in the formation of the 1985 DEM (Motyka et al., 2011).

Model domain
The forward model has to be evaluated repeatedly in the inversion, but all runs are instantaneous.This eliminates the need for a careful treatment of the boundary areas or the solution of the SSA in the entire drainage area, as done in regional time-dependent models.Instead we choose a limited model domain for efficiency, but include enough area around the used data sets (DEMs and bed elevation) to minimize boundary effects.We evaluate results spatially and along a centerline, which was extracted by approximately following the minimum bed elevation (Fig. 1). Figure 1 shows the model domain and the areas of high resolution DEMs and bed elevation as well as the misfit area used to calculate the datamodel misfit in the inversion.The SSA is solved over the entire model domain, but only velocity data within the misfit area is used to adjust the effective basal yield stress.Results are only interpreted within the misfit area, which is taken to be the same for all years.Areas outside the misfit area are shaded or excluded in all figures.

Choices in forward model and inversion
The model outlined above contains several poorly constrained parameter choices.In this section we discuss the choice for ice softness A in the forward model, the choice of model norm in the regularization term, the prior estimate for the effective basal yield stress, and the magnitude of the regularization parameter.For the model norm and the prior estimate of effective basal yield stress we used the 2006 data set, for all other parameters all inverted years where considered to determine the value.Final parameter choices were made after several iterations.We arrived at the following defaults values: Below we will discuss each choice by studying the effects of varying one parameter at a time, while holding the others at their default value.

Ice softness
The forward model contains many parameter choices, here we only discuss the ice-softness parameter.All other values for the forward model are discussed in Sect.2.1.1.Default values, or values that have proven to be good choices in other studies are used whenever possible.The SSA uses a viscosity that is dependent on a vertically averaged ice-softness parameter A which in turn depends on the temperature of the ice.Temperature has only been measured in a few boreholes (Lüthi et al., 2002) and its spatial distribution is not known.Here the vertically averaged ice softness does not vary horizontally for the entire model domain and we test different ice-softness values.A spatially variable ice softness would lead to effective basal yield stress fields that are consistent with the ice softness and therefore different than the effective basal yield stress fields found here.Nonetheless, we would expect all main findings about the changes and sensitivities of effective basal yield stress to stay true.Additionally, we conducted time-dependent numerical experiments (spin-ups), where not only the ice flow but also temperature fields were computed.These experiments show little horizontal variability in the vertically averaged ice softness.Suggested values of ice softness in Cuffey and Paterson (2010, chapter 3.4.6,p. 72 ff) range from 0.01 ×10 −24 Pa −3 s −1 for ice at −40 • C to 2.4 × 10 −24 Pa −3 s −1 for temperate ice, while values as high as 9.3 × 10 −24 Pa −3 s −1 have been reported from laboratory tests (Budd and Jacka, 1989).Higher values of ice softness are often used and justified by the anisotropy of ice or effects of grain size and/or impurities (Lüthi et al., 2002).
The achieved data-model misfit for different ice softness (Fig. 2) shows that only very hard ice (low A) leads to a marked increase in the data-model misfit.This confirms the finding of Joughin et al. (2012) that a hard ice model is not a good representation of the ice rheology of Jakobshavn Isbrae.On the other hand, Joughin et al. (2012) find with a terminus-driven model that a soft-ice model (A = 10 × 10 −24 Pa −3 s −1 ) does not transfer seasonal changes far enough inland.Here the ice-softness value 2.5 × 10 −24 Pa −3 s −1 is chosen for all years as a compromise between 2 and 3 × 10 −24 Pa −3 s −1 , which give the lowest datamodel misfit for 1985, 2000 and 2005, 2006, 2008.This ice softness is equivalent to an isothermal ice column with a temperature of ∼ −3 • C using the flow law temperature dependence given by (Cuffey and Paterson, 2010).For comparison, at a site on the ice sheet adjacent to the ice stream (Lüthi et al., 2002)  3) for different ice-softness values for all years.Hard ice (small value of ice softness A) leads to a marked increase in data-model misfit, whereas softer ice only slightly increases the data-model misfit.We choose the ice softness A = 2.5 × 10 −24 Pa −3 s −1 for all years and the range from 2-3 × 10 −24 Pa −3 s −1 is discussed.
ice, indicating our chosen ice softness has some enhancement relative to the borehole.

Model norm
The regularization term of the cost function contains a model norm (Eq.4).This term is necessary to stabilize the inversion.Choosing a model norm biases the solution and needs to be considered in the interpretation.As outlined in the Methods Sect.2.1.2,the type of model norm used here allows for a bias towards (1) "small" solutions, where the departure from a prior estimate of effective basal yield stress is penalized, (2) "smooth" solutions where the derivative of τ c −τ prior c is held small, which tends to preserve the shape of τ prior c , or (3) a mix between these two options.
Figure 3 shows L curves for three different model norms: pure L 2 norm, pure H 1 norm and L 2 norm with an additional small amount of H 1 norm.By increasing α more emphasis is placed on the data-model misfit minimization and more roughness is allowed in the solution.Calculating datamodel misfit values can be computationally expensive because each data point requires an inversion run and the inversions with very high α take many iterations to converge.
We show examples of modeled effective basal yield stress for under-and overfitting of the data, as well as a solution for the approximate "corner" of the L curve.The corner of the pure H 1 norm is at a data-model misfit approximately 50 m a −1 higher than the corner of the pure L 2 norm, and the effective basal yield stress field of the H 1 norm results in an accordingly smoother solution.but with similar limits for data-model misfits.This can be an indication of the total error, T tot , in the system.The pure L 2 norm produces large jumps in effective basal yield stress, especially with higher regularization parameter values, making it more sensitive to the choice of α.Here we prefer the pure H 1 norm solution because the non-localized nature of the SSA does not account for small-scale features in effective basal yield stress.Additionally, as long as the regularization parameter is chosen to yield similar data-model misfit values, the choice of norm influences the solution only within an acceptable range (see Sect. 5.1).

Prior estimate
In Tikhonov regularization, the cost function (Eq.2) is minimized, and a prior estimate of effective basal yield stress is necessary as a starting point for the iterations and for the model norm term.Within the misfit area the latter seems to outweigh the former.A prior estimate commonly used in glaciology is the driving stress field divided by two (Joughin et al., 2004).This choice was suggested because in the shallow-ice approximation the driving stress is locally balanced by the basal shear stress, but this is not necessarily the case for the SSA, where membrane stresses are considered.
Figure 4 shows two Tikhonov inversions with the prior estimate set to τ d /2 and to a constant value.Both of the resulting effective basal yield stress fields lead to almost identical residual velocity fields; in other words both solutions can account for the main features of the observed velocities.Small-scale features that are introduced in the τ d /2 prior estimate remain unchanged because they do not affect the ve- locity field sufficiently.The commonly used L 2 norm was applied for this figure; the H 1 norm would exacerbate the problem because the shape of the initial estimate tends to be preserved.
Without prior knowledge about the basal shear stress a constant prior estimate is most appropriate to avoid introducing small-scale features that may not be real.For the pure H 1 model norm adding a constant value to τ prior c will not influence the solution inside the misfit area.But we find that the inversion converges only for values within a certain range (approximately 5 × 10 4 -8 × 10 5 Pa).Therefore a good prior estimate could be the average of τ d inside the misfit area (here: τ d ≈ 1 × 10 5 Pa).Here we performed an inversion and used the value of modeled effective basal yield stress along the centerline at the upstream edge of the misfit area as the prior estimate (1.4 × 10 5 Pa).In this way the algorithm does not have to introduce extreme basal shear stress values to compensate for values outside the misfit area that lead to wrong ice velocities.All prior estimates in the remainder of this study were set to 1.4 × 10 5 Pa.

Regularization parameter
Given the choice of parameters discussed above, the L curve criterion can now be used to choose the appropriate regularization parameter, α.Commonly the L curve is displayed as a log-log plot, but for our inverse problem no clear corner emerges (not shown).There are different reasons for the lack of corner in the L curve, one of which is an increase in problem size (Hansen, 2001).As suggested by Calvetti et al. (2000) it is acceptable to plot data-model misfit against model norm on a linear scale to find the corner of the L curve. Figure 5 shows the linear plot of the L curve for all years with the above chosen parameters.We choose a regularization parameter of α = 10 for all years based on this figure.Data-model misfit values in Fig. 5 do not reach below 100 ma −1 , which is much higher than the expected root mean square error in surface velocity observations: assuming a 3 % error (Joughin et al., 2012) the root mean square error over the misfit area is ∼ 7 m a −1 .Errors are thus dominated by those introduced by the simplified model and/or geometry input fields.The high data-model misfit ensures that no overfitting of the observed surface velocity data occurs, but overfitting due to the model and parameter errors would still be a possibility without the regularization term.Since T obs is much smaller than the data-model misfit, we use the L curve method to improve parameters of the model such as the ice softness.

Results
Inversions for all years with the parameter choices discussed above are shown in Fig. 6.All inversions reproduce the overall pattern of observed surface velocities.This shows that, in general, our data and model choices are capable of reproducing the observations by only adjusting effective basal yield stress.But a small data-model misfit by itself does not speak to the quality of the resulting effective basal yield stress solution.
The first leg (lower 5 km of the glacier) shows a trend from higher to lower effective basal yield stresses over the years.Additionally, a slight widening of the area with low effective basal yield stresses is evident.The 2008 inversion results show continued widening, but the low effective basal yield stress area does not extend as far inland as in the 2006 results.Despite the use of independently produced DEMs and observed surface velocity data sets, the general spatial distribution of effective basal yield stress outside of the main fast flowing glacier remains fairly constant in all inversions compared to the large changes in the first leg.This consistency across years in areas with minimal observed changes in geometry and flow is encouraging and justifies the use of constant parameters for all inversion runs.Our main area of interest is the lower glacier with the largest changes in effective basal yield stress across the years.This area entails high values of observed surface velocities and a deep trough in the bed topography.Residual velocities (difference of modeled and observed) are generally high in this area of fast flow, but relative residuals are in fact similar or lower than in the slow flowing areas (Fig. 6).
To compare the results for the different years in more detail, Fig. 8 shows the results along the centerline for all years.Here the basal shear stress, τ b , calculated according to Eq. ( 1), is shown and compared to the driving stress.As seen in the spatial distribution of effective basal yield stress, the values in the first leg are clearly lowered compared to higher upstream, and they generally decrease over time.Despite minimal changes in driving stress from 1985 to 2006, the basal shear stress changes significantly over this time period.In 2000 only a lowering close to the first bend is visible, whereas basal shear stress close to the terminus increases compared to 1985.Past the first bend, the inverted basal shear stresses are generally higher; for 2008 the average value of τ b in the first leg is 0.2 × 10 5 Pa, whereas the average value between the first and the second bend is 1.8 × 10 5 Pa.Upstream of the first leg no clear trend in basal shear stress is visible, which is in contrast to the general increase in basal shear stress in this area inferred by Joughin et al. (2012).The basal shear stress accounts for about 20-40 % of the driving stresses along the entire centerline, with a few single peaks reaching 80-100 % of the driving stresses.

Robustness of inversion
The solution to our inverse problem is not unique, many of the parameters are not well constrained and a range of parameter choices would be equally acceptable.The emphasis here is on temporal changes in effective basal yield stress, and little significance should be given to the actual value of the stress in a given inversion.To evaluate the robustness of our results, we explore a range of parameters for the years 1985 and 2006.
We chose an ice-softness value of 2.5 × 10 −24 Pa −3 s −1 for all years, while the minimum data-model misfit values are reached for ice-softness values between 2 and 3 × 10 −24 Pa −3 s −1 (see Fig. 2). Figure 9 shows an envelope of solutions of effective basal yield stress along the centerline for this range of ice softness.The solutions for A = 3×10 −24 Pa −3 s −1 lead to generally higher τ c values than the A = 2×10 −24 Pa −3 s −1 solutions, because softer ice leads to a more localized stress balance and therefore to higher values in effective basal yield stress.The 2006 effective basal yield stress solution exhibits a higher sensitivity to changes in ice softness and the effective basal yield stress is affected most just upstream of the first bend.It is important to keep in mind that we are using a constant value of ice softness over the entire model domain.Larger variations of effective basal yield stress are possible for more realistic representations of the temperature distribution in the ice.As a thermomechanically coupled ice-sheet model, PISM is capable of producing realistic ice temperature fields, which could be achieved through spin-ups.But it is not clear which effective basal yield stress values to use for such a spin-up.Joughin et al. (2009) for example used iterative spin-ups to find an ice temperature field that is consistent with the effective basal yield stress.
One of the most important sources of uncertainty is the choice of regularization parameter.As mentioned before, it is not straight forward to choose the exact location of the "corner" in the L curve.In other studies the regularization parameter is chosen by calculating the point of maximum curvature (Vogel, 1987, ch. 7.4).But even when this point is calculated exactly, the L curve criterion remains an approximate method.Therefore, we chose the approximate value of α = 10 and an upper and lower bound (α = 3 and α = 30).Figure 10 shows that the choice of regularization parameter mostly affects the first leg where a smaller data-model misfit in velocities is expensive (in the model norm sense) because the narrow trough makes abrupt changes in τ c necessary.The data-model misfit is a root mean square over the misfit area, meaning that local under-or overfitting is possible (and very probable).When plotting the data-model misfit relative to u obs along the centerline for different regularization parameters (Fig. 10), it becomes clear that in the first leg the fit to velocity observations is still improving, unlike in areas higher upstream.A higher α could be justified when focusing on the inversion results of τ c in the first leg.
We also want to investigate how a different choice of model norm would have affected our solution.For a direct comparison with the range of regularization parameters used for Fig. 10, we chose a conservative α = 0.01 as the "ideal" solution and a range from α = 0.003 to α = 0.03 (Fig. 11).The sharp features in τ c for α = 0.03 between the two bends reach values of 4.1×10 5 and 2.6×10 5 Pa for 1985 and 2006, respectively, showing the sensitivity of this norm to overfitting.Note that the relative residual does not improve significantly even though such large features are introduced.The actual corner of the L 2 L curve is at α = 0.1 and the solution for this regularization parameter is shown as well in Fig. 11.The value of the modeled τ c is generally lower in this case and displays sharper features, while the improvement in relative residual is not significant.
To illustrate how a prior estimate with small-scale features can influence the solution, Fig. 12 shows the centerline solutions for a prior estimate of τ d /2.When using prior estimates with small-scale features, the L 2 norm is more useful because it does not try to conserve the shape of the prior estimate.The centerline solution only contains fast flow, where τ c is adjusted well, in slow flow areas there are more places where the small-scale features of the prior estimate remain.Half of the driving stress might be a good first order approximation of effective basal yield stress, but when applied unsmoothed as a prior estimate, it introduces spurious features.To initialize entire or drainage-basin-wide ice-sheet models a continuous field of effective basal yield stress is needed.The inversion algorithm only calculates the data-model misfit where surface velocity observations are available.This can lead to large areas where the prior estimate will determine the final τ c .Future work should consider what the best strategies for the prior estimate in such situations are.To compare inverse results of different years, ideally we would use inverse methods where the cost function also includes a penalization for changes in time as done in time-dependent seismic tomography (Julian and Foulger, 2010).In this manner years with larger velocity data coverage would adjust τ c in areas with data gaps in other years.

Changes in effective basal yield stress
Figure 6 shows a general decrease in effective basal yield stress close to the grounding line, here we explore how this relates to changes in geometry.We solely concentrate on snapshots of ice geometry and do not investigate causes of the change in geometry, such as increased melt or decreased buttressing at the ice front.In other words, the inversion examines an instantaneous stress state given a certain geometry and surface velocity, but it cannot, by itself, attribute any causes.A common way to parameterize the effective basal yield stress in time dependent model runs is through a Mohr-Coulomb model (Iverson et al., 1998): where (ρgH − p w ) is the effective pressure, p w the pore water pressure, g the gravitational acceleration, ρ the density of ice (set to 917 kg m −3 ), and φ a "till friction angle," a strength parameter for the till comparable to "angle of repose" for granular piles.To find out if the changes in the inverted τ c are in agreement with such a parametrization, we compare the relative change in τ c (LHS of Eq. 6) to the relative change in height above floatation (RHS of Eq. 6).We assume that the basal water pressure is equivalent to oceanic pressure (p w = ρ w g|z b |, ρ w is the density of water and set to 1025 kg m −3 , where the bed elevation is below sea level and p w = 0 otherwise).The term tan(φ) cancels when calculating the relative change (e.g., for 1985 and 2006): The proximity to floatation is important in this calculation and we are subtracting 30 m (approximate offset at the 2007 terminus) from |z b | to correct for the geoid-ellipsoid separation in the area of terminus.The area of interest lies entirely in the ablation area, so that density variations due to firn do not need to be considered.Density variations caused by heavy crevassing, however, can occur, but are not considered here.
Figure 13 shows that the relative change in inferred τ c is much more localized to the trough than the relative change in height above floatation.A slight increase in τ c is visible near the margins of fast flow.But the broad pattern is similar, confirming that the relative change in height above floatation accounts for most of the relative changes in τ c .Also for shorter timescales and after the disintegration of the floating ice tongue (2005)(2006) similar patterns of relative change are visible (Fig. 13).An increase in sliding due to more melt water at the base, for example, would lead to a spatial pattern of relative change distributed over the entire area of melt.Because we do not see this spatial pattern related to melt area, our results support the findings of Joughin et al. (2008a) that increase in seasonal melt is not the main driver of the observed speed-up.
Figure 14 shows how the relative change in inferred τ c (LHS of Eq. 6) and the predicted relative change in height above floatation (RHS of Eq. 6) compare along the centerline.The relative changes in inferred τ c are shown for a range of regularization parameters and ice softness.The relative change in height above floatation has a different qualitative shape, but falls within the envelope of regularization parameters.The choice of regularization parameter gives a large uncertainty in relative changes in τ c , especially in the terminus area.Above we showed that there is a significant lowering in τ c in the first leg, even when taking into account the uncertainties introduced by the parameter choices in the inversion.Figure 14 on the other hand, shows that these same uncertainties of the inversion method make it difficult to judge the validity of parameterizations for τ c (Eq. 5).
To investigate if using a constant-in-time value for the till friction angle φ is reasonable, we plot the inferred value of τ c against the predicted effective pressure for each grid point.In areas with a constant till friction angle we would expect a linear relationship between τ c and effective pressure (ρgH − p w ) with a slope of tan(φ).The overall thinning from 1985 to 2006 should lead to a decrease in effective pressure and a simultaneous decrease in τ c .We expect the same linear relationship for both years, but with a data point cloud shifted towards lower values of effective pressure for   2006.When taking into account the entire misfit area, no relationship is apparent (Fig. 15), but when we limit the analyzed points to the areas of fast flow, a linear relationship emerges (Fig. 15).The slope of this linear fit indicates that tan(φ) ≈ 0.02 and thus φ ≈ 2 • , which is a very low value of till friction angle compared to the measured values between 19 • and 26 • (Iverson et al., 1998;Kamb, 1991).The consistent linear relationship in the fast flow area and the shift in data points to lower values are in agreement with the assumption of a constant tan(φ) in time.The unphysical value of φ and the lack of relationship between τ c and the effective pressure over larger spatial scales, however, show that a simple parameterization might not adequately represent the actual bed properties under Jakobshavn Isbrae.In this study we use an approximation to a perfectly plastic sliding law, there- fore, τ c is only an approximation to the basal yield stress.We test additional smaller values of q (q = 0.1 and q = 0.001) in Eq. ( 1) to see if a closer approximation to a plastic till affects our findings.The actual value of τ c increases by up to 2 × 10 5 Pa, but the lowering of τ c in the first 7 km during the time of acceleration is a robust result.Comparing τ c to effective pressure leads to slightly higher values of till friction angle (φ ∼ 3 • ), but these values are still low compared to measured values mentioned above.The inversion calculates the best fit to observed velocities, given an SSA forward model and the restrictions from the regularization.If results from the inversion are used in prognostic forward models that are based on the SSA, it might be more appropriate to use these inversion values, even if they differed significantly from actual in situ measurements of till friction angle or effective basal yield stress (if those were indeed measurable).In that sense, the goal of an inversion is not always to find the true physical parameters, but rather those that are consistent with a simplified physical model and the observations.

Conclusions
A careful choice of parameters in an inversion is especially important when comparing effective basal yield stress distributions independently inferred for different years.To estimate the influence of the parameter choices, reasonable ranges are explored and we find that the weakening of effective basal yield stress over the years close to the terminus area is a real temporal variation.The observed changes are in agreement with a Mohr-Coulomb parameterization of effective basal yield stress, where the change in effective pressure is the main driver for the changes in effective basal yield stress.Despite this broad agreement, the involvement of other processes cannot be excluded and the sensitivity of the inversion to parameter choices, in particular the regularization parameter, makes it difficult to evaluate effective basal yield stress parameterizations.The spatial distribution of residuals shows that for Jakobshavn Isbrae less simplified models, improved bed topography and/or a spatially varying ice softness could potentially improve the inversion results.With the currently available satellite data and the length of observational record on many other fast changing glacier systems it is possible to apply these methods to other systems and to further advance our understanding of the changes at the base of the ice.

Fig. 1 .
Fig. 1.Model domain (entire area shown) with MODIS image for reference (single pass MODIS image, spring 2001, courtesy of M. Fahnestock).Also shown are the extent of the higher resolution bed topography (cyan), 2007 DEM (green), 1985 DEM (blue), misfit area (red), straightened centerline (dashed black), the white circles mark the two "bends" that are mentioned later on and the area later shown in all the map-view figures (solid black).

Fig. 2 .
Fig.2.Data-model misfit M (Eq. 3) for different ice-softness values for all years.Hard ice (small value of ice softness A) leads to a marked increase in data-model misfit, whereas softer ice only slightly increases the data-model misfit.We choose the ice softness A = 2.5 × 10 −24 Pa −3 s −1 for all years and the range from 2-3 × 10 −24 Pa −3 s −1 is discussed.

Fig. 4 .
Fig. 4. Influence of different prior estimates.Inversions with 2006 velocity data and a prior estimate of effective basal yield stress of (Top) τ prior c = τ d /2 and (Bottom) τ prior c = 1.4 × 10 5 Pa.The columns show the prior estimate, the inferred effective basal yield stress, the change of the prior to the modeled τ c and the residual in velocity (|u obs − u mod |).Pure L 2 model norm, α = 0.1.

Fig. 5 .
Fig. 5. L curves for all years plotted on a linear scale.The range of regularization parameters is α = 0.1 -1 × 10 3 .Based on this figure α = 10 is chosen for all years.

Fig. 6 .Fig. 7 .
Fig. 6.Inversion results for 1985, 2000, 2005, 2006 and 2008.The columns show the modeled τ c (logarithmic scale), the velocity residual (|u obs − u mod |), the relative velocity residual (100 |u obs − u mod|/u obs ), the observed velocities u obs and the modeled velocities u mod .The area past the 2008 grounding line is not included in the misfit area and is blacked out.

Fig. 8 .
Fig. 8. Inferred basal shear stress, τ b , along centerline for all years.Area outside of misfit area is shaded gray and the blue vertical lines show the position of the two bends in the centerline.(Top) Crosses mark the driving stresses, τ d , for the years 1985 and 2006.The sharp peak in τ b occurs at the grounding line for each year.(Bottom) Modeled (solid lines) and observed (points) velocities for all years.

Fig. 11 .Fig. 12 .
Fig. 11.Robustness of effective basal yield stress results for L 2 norm with the conservative regularization parameter values α = 0.003 (upper envelope) and α = 0.03 (lower envelope), the black line indicates the α = 0.01 solution.The actual corner of the L 2 L curve is at α = 0.1 and the solution for this regularization parameter is shown in cyan.

Fig. 13 .
Fig. 13.Relative change in inferred τ c (left) compared to the change predicted by the Mohr-Coulomb parameterization used in PISM.Areas where the bed topography, z b , is above sea level are masked out.

Fig. 14 .
Fig. 14.Relative change in inferred τ c along the centerline for a range of regularization parameters (top) and ice softness (bottom).The red line shows the change in height above floatation predicted by the Mohr-Coulomb parameterization used in PISM.(a) compares years 1985 and 2006, (b) compares the years 2005 and 2006.

Fig. 15 .
Fig. 15.Inferred τ c against effective pressure (ρgH − p w ) for each grid point where the observed velocities are greater than the threshold velocity given above the inset plots.A linear fit is given for both years, the slope is 0.016 (0.019) and the intercept is 1.37 × 10 4 Pa (0.97 × 10 4 Pa) for 1985 (2006).