Interactive comment on “ An iterative inverse method to estimate basal topography and initialize ice flow models ”

This paper presents a method for recreating the bedrock topography of a glacier given that (a) the surface topography is known for a specific point in time, (b) the mass balance history can be estimated for a sufficiently long time frame, and (c) a model for the ice dynamics is available. The method consists in running the ice-dynamics model forward in time, and iteratively adapting an initial guess of the bedrock in order to minimize the mismatch between modeled and known surface geometry. The merit of the paper doesn’t lie much in the application in three dimensions of a methodology which is already found in the literature (see Heining, Phys. Of Fluids, 2011, or Michel et al, Inv Proc., 2013), but in the clarity of the presentation, making the method very easy to C222


Introduction
The inaccessibility of the bed underneath glaciers and ice caps limits observational constraints on basal conditions.Surface data are much more abundant and, when linked to basal parameters using a forward model, can be used in an inverse manner to better constrain basal boundary conditions (e.g.MacAyeal, 1992;Joughin et al., 2004;Raymond Pralong and Gudmundsson, 2011).In this study, we focus on reconstructing basal topography.Accurate estimates of the basal topography and ice thickness distribution are needed mainly for two reasons.Firstly, estimates of the ice thickness distribution are required to determine the total ice volume stored in a glacier or ice cap, which is needed to assess the potential contribution to sea-level rise.Secondly, accurate simulation of the future evolution of glaciers and ice caps in a changing climate depends strongly on the initial ice thickness distribution.Proper initialization of ice flow models should result in an ice geometry and thermodynamic state at the start of a prognostic experiment, which is consistent with the history of the surface climate forcing.The relevance of accurate initialization has been stressed by, for example, Oerlemans et al. (1998), Huybrechts and de Wolde (1999) and Huss et al. (2008).
Measuring distributed fields of basal topography using ground penetrating radar (GPR) is a laborious and expensive task, hence bed height data sets generally cover the region of interest only partially.This has motivated the development of techniques to indirectly estimate the ice thickness distribution.A first attempt was made by Nye (1951), who assumed a perfectly plastic ice rheology in order to compute the ice thickness from the surface slope and a constant yield stress.This method has been used to recover glacier bed topographies by, for example, Haeberli and Hölzle (1995), Oerlemans (1997) and Paul and Svoboda (2010).Recently, Li et al. (2012) extended this method accounting for variable Published by Copernicus Publications on behalf of the European Geosciences Union.
width of the glacier, thereby incorporating the effect of side drag on the stress state.Alternatively, several studies used the mass continuity equation to derive the thickness distribution.Farinotti et al. (2009) and Huss and Farinotti (2012) applied a mass conservation approach to compute the thickness distribution from surface height data.Other mass continuity conservation approaches rely on velocity data to compute bed heights (e.g.Morlighem et al., 2011;McNabb et al., 2012).Inverse control methods have been applied to Antarctic ice streams to simultaneously recover basal topography and slipperiness, requiring distributed fields of surface topography and velocity (Thorsteinsson et al., 2003;Raymond and Gudmundsson, 2009;Raymond Pralong and Gudmundsson, 2011).In Clarke et al. (2009), an alternative method is proposed and discussed, which uses neural networks to obtain thickness estimates from surface height data.
In this study, an iterative inverse method is presented to reconstruct distributed basal topography, given a map of surface height data, an ice flow model and a time-dependent surface climate forcing.The approach does not require surface velocity input and aims to minimize the discrepancy between modeled and observed surface heights at the end of timedependent model runs.In this way, it additionally provides an efficient way to initialize ice flow models for forecasting experiments.In the spirit of a recent studies by Pollard and DeConto (2012) and Michel et al. (2013), a simple regularization method is employed, in which the surface height misfit is directly used to iteratively update basal topography.In this manner, complex evaluation of cost-function gradients in order to define a search direction is avoided.This study builds on previous applications of such a correction method in a flow-line context (Oerlemans, 2001;Leclercq et al., 2012;Michel et al., 2013) and extends the analysis to three dimensions.The real glacier application to Nordenskiöldbreen, with added complexities related to the presence of data errors and climate variability/feedbacks, serves as a step-by-step manual for future applications recovering distributed basal topography and/or initializing ice flow models.First, we present the regularization method and discuss relevant features of the forward model (Sect.2).Then the approach is applied in a set of synthetic steady-state experiments to evaluate convergence and robustness of the inverse method in three dimensions (Sect.3).In Sect. 4 the inverse approach is applied to the glacier geometry of Nordenskiöldbreen, Svalbard.Reconstructed beds are validated against radar observations and the significance of errors is discussed.

Model description
The Parallel Ice Sheet Model (PISM) is used as a forward model in the iterative procedure in which the basal topog-raphy is adjusted after every iteration.Here we will address relevant features of the model.For a more elaborate description of the model, the reader is referred to Bueler and Brown (2009) and Winkelmann et al. (2011).
Forced with fields of surface mass balance and nearsurface temperature, experiments are performed with PISM to simulate the spatiotemporal evolution of the surface topography.Ice flow in PISM is a combination of ice movement by internal deformation, described by the shallow ice approximation (SIA) (Morland and Johnson, 1980;Hutter, 1983), and basal sliding, described by the shallow shelf approximation (SSA) (MacAyeal, 1989;Weis et al., 1999).Weighting and averaging of the two velocity components results in a smooth transition from shearing to sliding flow and hence provides a solution to the flux-matching problem at the transition from sliding to non-sliding ice (Schoof, 2006;Bueler and Brown, 2009).Effectively, the SSA is used as a sliding law.As a direct consequence of its shallowness, the model is significantly more efficient than higher-order or full-Stokes alternatives in terms of computational costs.
The evolution of the internal thermodynamic structure of the ice is modeled using an energy conservation scheme in which an "enthalpy-gradient approach" is adopted to compute vertical profiles of water content and temperature (Aschwanden and Blatter, 2009;Aschwanden et al., 2012).Enthalpy advection, diffusion and production by strain heating are modeled in a consistent manner in both cold ice (temperature below pressure-melting point) and temperate ice (temperature at pressure-melting point).Additionally, a bedrock model simulates the diffusive enthalpy transport in the upper hundreds of meters of bedrock underlying the ice.A fixed geothermal heat flux (0.042 W m −2 ) is prescribed at the lower boundary of the bedrock model.When ice slides over the bed, frictional heating will contribute to the heat budget at the ice-bed interface.In case the bed is at pressuremelting point, excess energy will induce basal melting and local water production.Basal water is stored locally and may refreeze when the basal temperature drops below pressuremelting point.In this simplified water model exchange of water between grid cells is not accounted for.An approximation of the local conservation of energy determines the spatial distribution of temperate and cold ice and the spatial distribution of basal water production.A more realistic water model, incorporating a Darcian transport scheme, surface to bed transport and a transient transmissivity of the system, would be needed to realistically simulate water transport and seasonal variability in water layer thickness and pressures, which is beyond the scope of this work.The connection between water layer thickness, water pressure and the basal shear stress follows the description of Bueler and Brown (2009) and Van Pelt and Oerlemans (2012).A pseudo-plastic formulation for the basal shear stress τ b as a function of the basal sliding velocity u b and the yield stress τ c is used: A value of the pseudo-plasticity exponent q of 0.4 is used here, which is equivalent to a Weertman-type sliding law exponent of 2.5.The yield stress τ c depends on the material till strength φ and water pressure p w according to Clarke (2005): where ρ is the ice density (900 kg m −3 ), g the gravitational acceleration (9.81 m s −2 ) and H the ice thickness.The water pressure p w is a fraction of the overburden pressure and depends linearly on the effective water layer thickness as described in Bueler and Brown (2009).Given the simplicity of the determination of the basal shear stress, which is inevitable given our current knowledge of the basal conditions, we accept that a major source of error in the reconstructed basal topography stems from uncertainty in modeling the basal resistance.In Sect. 4 reconstructed beds will be validated against independent ground-penetrating radar data, which enables us to quantify the accuracy of the reconstruction and assess the significance of modeling errors.

Correction method
Our main aim is to develop, test and apply an inverse approach in which bedrock heights are iteratively adjusted in such a manner that the misfit of modeled and observed surface topography is minimized.Here we discuss the regularization technique used to stabilize the inverse procedure.In contrast to many previous inverse approaches in glaciology, often aiming to minimize the misfit in surface velocities to reconstruct basal shearing conditions (e.g.(MacAyeal, 1992), (Joughin et al., 2004)), the correction method used here does not rely on evaluation of cost function gradients to update model parameters.Pollard and DeConto (2012) recently presented a correction method in which the local surface height misfit is directly used to iteratively update basal slipperiness.
In a similar fashion, we directly use the surface height misfit to iteratively adjust bed heights.Forced at the surface with prescribed mass balance and near-surface temperature input, the ice flow model (PISM) is iteratively run over a specified period of time.At the end of every model iteration, local deviations between modeled surface heights h mod and reference heights h ref are evaluated and directly applied to compute a new bed b: where K is a constant, referred to as the relaxation parameter, and n is the iteration number.The magnitude of the bed correction in response to a surface discrepancy hence scales linearly with the relaxation factor K, which needs to be chosen small enough to avoid instabilities due to overcompensation of the bed, and large enough to speed-up convergence of the approach.The same correction method has previously been applied in a flow-line context in Oerlemans (2001), Leclercq et al. (2012) and Michel et al. (2013) and is based on the intuitive principle that a higher bedrock elevation leads to an overall higher surface elevation.In the field of fluid dynamics Heining (2011) studied gravity-driven film flow over an undulating bed and used the correction method to successfully recover basal topography in idealized flow-line experiments.Perturbations of the bed height lead to a non-local surface expression with a transfer amplitude and phase-shift depending on multiple factors including slip ratio and dimensions of the basal perturbation, as shown by Gudmundsson (2003) and Raymond and Gudmundsson (2005).Observational verification for many of their theoretical findings has recently been provided by De Rydt et al. (2013).When recovering a single bump in an otherwise flat bed, the surface misfit and hence the bedrock correction will initially be out-of-phase with the actual bump position (when starting from a fully flat bed).We will illustrate in Sect.3.1 that the interaction of bed adjustments and the surface misfit causes the recovered bump to gradually grow and move into the right position.

Synthetic experiments
Experiments with a synthetic ice sheet give insight in the functionality of the bedrock reconstruction method and the sensitivity of convergence to the chosen setup of the experiments.
Prior to the iterative procedure to recover basal topography, reference bedrock and surface height geometries are required.The reference geometries are the geometries we aim to reconstruct with the inverse procedure.A reference bed topography is manually designed and serves as a lower boundary in an initial model run performed to produce a reference surface (h ref in Eq. 3).After the initial run, the designed reference basal topography is "forgotten" and the iterative recovery starts from a random initial bed (a flat bed in the synthetic experiments) and attempts to gradually recover the reference basal topography.In both the initial and iterative model runs, the model starts without ice and an invariant mass balance and temperature forcing is applied, shown in Fig. 1.The model runs forward in time over a period of 20 kyr, which is sufficient to reach a steady-state geometry.After every model iteration, the surface height misfit (h mod − h ref ) is computed and bed heights are adjusted accordingly (Eq.3).Both the prescribed mass balance and temperature at the snow-ice interface are constant in time and vary linearly with distance to the domain center (Fig. 1).Resulting steady-state geometries are approximately axisymmetric and the velocity field is characterized by an internal deformation dominated center region and more significant sliding velocities towards the ice margin.In the standard setup the relaxation factor K is set to 1.
Experiments are conducted on a 50 × 50 horizontal grid with a 15 km grid spacing.The vertical grid in the ice contains 100 unevenly distributed layers (more concentrated near the bed).The vertical grid in the bedrock consists of 10 layers extending down to 200 m below the ice-bedrock interface.
The idealized framework of the synthetic experiments allows us to study convergence of the inverse approach in three dimensions under simplified conditions.Additional complexities related to the presence of noise in the surface height data, height feedbacks of the climate forcing and nonsteadiness of the final surface profile, which are relevant in an application with real data, are addressed in Sect. 4. The experiments in this section are divided into two parts, discussing recovery of a bed containing a single bump ("onebump" experiment, Sect.3.1) and a distributed field of bumps ("bump-field" experiment, Sect.3.2).

"One-bump" experiment
In a first experiment the iterative procedure is used to recover a single bump in an otherwise flat bed.The initial run produces reference surface height and bedrock height profiles shown in Fig. 2a and b, respectively.In the iterative bed recovery procedure, an arbitrary total of 40 iterations is performed of which snapshots of the recovered bed are shown in Fig. 2d.Corresponding snapshots of surface misfit, defined as the modeled minus the reference bed height, are plotted in Fig. 2c.Surface velocities are dominated by deformation in the interior (up to 13 m a −1 ) and sliding near the ice margin (up to 48 m a −1 ).
With the increasing number of iterations the bump in the bed is gradually recovered (Fig. 2d) and the surface misfit diminishes (Fig. 2c).The presence of the bedrock bump in the reference bed raises the reference surface upstream of the bump position (and slightly lowers the surface downstream of the bump).At the start of the iterative procedure with a flat bed, the surface misfit is therefore negative upstream of the reference bump position, thereby inducing a positive bed correction.In subsequent iterations the gradually growing bump upstream of the reference bump position affects the surface misfit such that the negative misfit (positive bed correction) gradually shifts in the downstream direction (Fig. 2c).Consequently the recovered bump slowly moves towards the reference position, while its amplitude increases (Fig. 2d).It can be shown that for a phase shift between the reference bed bump and the corresponding surface bump smaller than 180 • , the recovered bump position will converge to the reference position.Raymond and Gudmundsson (2005) and Gudmundsson (2003) illustrate that phase shifts do not become larger than 90 • for a wide range of parameter values.A crosssectional view of the bed in Fig. 2e further illustrates the recovered bump movement and growth with the number of iterations.The above experiment demonstrates that despite the non-local surface expression related to a local bed height perturbation, the local recovered bed approaches reference values when applying sufficient iterative corrections.The nonlocal surface expression in both the along-flow and crossflow directions is apparent for the entire grid (Fig. 2c) and reduces with the number of iterations.Note that the bed perturbation of up to 150 m results in a surface perturbation of only up to 10 m (ice thickness ∼ 1000 m), which is indicative of the low transfer amplitude of bed height variability to the surface.A larger transfer amplitude would enhance the growth rate of the recovered bump and hence speed-up convergence.Among other factors, the low mean surface slope (0.3 • ) in this example limits the bed to surface transfer amplitude (Raymond and Gudmundsson, 2005).

"Bump-field" experiment
The previous experiment demonstrates convergence of the bedrock reconstruction method when recovering a single bump.Next, we apply the iterative method to a distributed field of bumps and additionally verify the robustness of the approach when changing (1) the initial bed, (2) the relaxation factor K, (3) the mass turnover and (4) the bump dimensions.
An initial model run produces reference bed and surface heights shown in Fig. 3a, and b.Snapshots of the bed during the iterative procedure in Fig. 3c indicate gradual reduction of the misfit of the reconstructed bed with regard to the reference bed.A cross-sectional profile of the bed in Fig. 3d illustrates convergence with the number of iterations and shows that after 40 iterations, errors remain largest near the center of the ice sheet.Raymond and Gudmundsson (2005) showed that the bed to surface transfer amplitude after a bed height perturbation is proportional to the surface slope and the slip ratio, i.e., the ratio of sliding vs. shearing velocity.In this experiment, enhanced sliding velocities towards the margin and low slopes in the interior lead to an increasing transfer amplitude towards the margin.A higher transfer amplitude implies a stronger surface response to bed height deviations and consequently larger corrections of the bed.From this, lower bed height convergence rates towards the ice sheet center in Fig. 3c and d can be understood.
Figure 4 shows the root-mean-square surface misfit versus the number of iterations for multiple sensitivity experiments.The magnitude of the bed correction scales linearly with the relaxation factor K. A low value for K implies slower convergence of the bed.On the other hand, setting K too high might lead to overcorrection of the bed.That is, when the surface height is very sensitive to changes of the bed height, the associated large bed corrections might lead to a surface misfit (of opposite sign) larger than its value after the previous iteration, which is a sign of overcompensation of the bed.Consequently, an optimum value for K can be selected, which is small enough to achieve stable convergence and large enough to increase the rate of convergence.Fig. 4 clearly shows faster convergence for large K.
The sensitivity of the approach to a different bed topography at the start of the first iteration is investigated by performing an additional run in which the initial bed height is strongly overestimated by setting it to the reference surface height.This will lead to a significant overestimation of the surface height in the first few iterations.Fig. 4 illustrates that, despite a much larger initial bed (and surface) misfit, the misfit after 40 iterations is nearly identical.Final beds after 40 iterations are similar, indicating that the reconstructed bed is not very sensitive to the chosen initial topography after many iterations.
We further investigate the applicability of the method when applying a different mass balance-altitude gradient, implying significant changes in ice dynamics and surface geometry.Before starting the iterative recovery, reference profiles are constructed consistent with the new surface forcing.Figure 4 illustrates that convergence is apparent for both a  positively (4 times larger) and negatively (4 times smaller) perturbed balance gradient.A larger balance gradient implies larger ice fluxes and a thicker ice mass with steeper slopes near the margin.Steep slopes in combination with a high slip ratio near the margin result in a high bed to surface transfer amplitude.The wavy pattern of the surface misfit for the large balance gradient can be ascribed to overcorrection of the bed near the margin.Setting K to a smaller value would reduce this effect.
In order to study the sensitivity of convergence for varying bump dimensions, the inverse approach is applied to recover a bed containing a mixture of two bump sizes (Fig. 5b).Snapshots of the bed in Fig. 5c and the bed height misfit series in Fig. 5d indicate convergence and illustrate that large-scale features are recovered more quickly than smaller-scale features.The recovery of the large-scale bumps leads to a rapid drop in the bed misfit until around 15 iterations.The slower recovery of the small-scale bumps continues until the end of the experiment and is particularly slow in the interior (again as a consequence of the small slope).The slower recovery of smaller sized bumps is in line with studies by Gudmundsson (2003) and Raymond and Gudmundsson (2005), which show an inverse dependence of bed to surface transfer am-plitude with the frequency of the bed perturbation.The fact that larger features are recovered prior to smaller features is used in Sect. 4 to define a stopping criterion for the iterative procedure, which is required when dealing with actual data (containing noise) and modeling errors.
In the aforementioned experiments reference surface heights were generated with the model prior to the inverse procedure.In that case a bed is known to exist for which the reference surface heights exactly match modeled values.Therefore, the recovered bed continues to improve with the number of iterations.Nevertheless, we hypothesize there is a numerical limitation to the detail in the bed that can be recovered because beyond some point very small bed adjustments may no longer affect surface heights due to numerical rounding/diffusion.In real world applications, surface height data contain a certain amount of error and the model forcing and physics are imperfect.In that case a perfect match between modeled and observed surface heights is not desired.When the iterative procedure is continued for too many iterations, further reduction of small-scale surface errors induce unrealistic adjustments of the recovered bed.Beyond a certain amount of iterations, adjustments of the bed based on the remaining surface misfit are hence no longer an improvement.As discussed in the next section, early termination of the inverse procedure can be applied to prevent contamination of the bed with undesired noise.

Application to Nordenskiöldbreen, Svalbard
The experiments in the previous section have illustrated that the inverse approach enables recovering basal topography in an idealized setting.Next, the applicability of the approach in a real world situation is investigated.The inverse method is applied to the geometry of Nordenskiöldbreen, Svalbard, for which we have an accurate set of surface height data from a digital elevation model (DEM) (Korona et al., 2009) and some bed height data, gathered using GPR.The DEM serves as a reference surface, whereas the bed height data are used for validation of the reconstructed bed.Compared to the idealized examples in the previous section application to Nordenskiöldbreen poses some additional challenges.First of all, as mentioned before, the presence of data and modeling errors suggest the need for termination of the iterative procedure before fitting the modeled surface to unrealistic features.Hence, a stopping criterion is used to minimize the significance of overfitting.Secondly, in contrast to the synthetic experiments, the glacier is not in steady state at the end of an iteration.A time-dependent surface climate forcing since 1300 AD is prescribed and the feedback effect of surface height changes on the mass balance and surface temperature is taken into account, as outlined in Sect.4.2.

Nordenskiöldbreen
Nordenskiöldbreen is a tidewater glacier in central Svalbard, connected to a large ice plateau, Lomonosovfonna.After recent retreat the ice front has partly retreated on land and the significance of frontal calving has declined (Plassen et al., 2004).GPS observations of ice movement along the flow line since 2006 reveal annual mean surface velocities up to 60 m yr −1 (Den Ouden et al., 2010).Nordenskiöldbreen has a polythermal thermodynamic structure with cold annual mean near-surface temperatures in the ablation area and a nearly temperate snow pack in the accumulation zone as a result of significant melt water refreezing (Van Pelt et al., 2012).This type of structure is common for Svalbard glaciers (Pettersson, 2004).A 40 m resolution DEM of the glacier and its surroundings, obtained in 2007 as part of the SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: Reference Images and Topographies) project (Korona et al., 2009), is shown in Fig. 6a.Some post-processing has been done to remove inconsistent data.The black line in Fig. 6a marks the extent of the selected grid of Nordenskiöldbreen, including two side-flows between Terrierfjellet, Ferrierfjellet and Minkinfjellet, merging with the main flow further downstream.GPR bed height data along the track shown in Fig. 6b, gathered using a 500 MHz radar during a field campaign in 2010, reveal bed heights below sea level upstream of the current frontal position.The GPR data are used to validate the reconstructed bed.

Climate forcing
In this section, the strategy employed to obtain a timeand altitude-dependent surface forcing is described.Given a winter air temperature reconstruction for Svalbard Airport, Longyearbyen (Divine et al., 2011), and an annual accumulation record for Lomonosovfonna, a climate reconstruction for Nordenskiöldbreen back to 1300 AD has been made.The climate sensitivity of the mass balance and subsurface temperature, constituting the surface forcing, is computed with a coupled energy balance-snow model.
The winter air temperature record is used in combination with temperature data from the Longyearbyen weather station to obtain annual temperature estimates.The instrumental record for Longyearbyen since 1912 reveals a much stronger year-to-year variability in winter than in summer temperatures.Additionally, over the instrumental period summer temperatures show a substantially smaller positive trend (0.011 K yr −1 ) than winter temperatures (0.024 K yr −1 ).In line with this, future climate scenarios for this region predict a higher sensitivity of winter temperatures in a changing climate (Førland et al., 2009).We compute the mean ratio of trends in winter temperatures and annual mean temperatures over the instrumental period (since 1912) and subsequently use this ratio to convert the winter temperature record, derived from a composite δ 18 O record from two ice cores in Svalbard (Divine et al., 2011), into annual mean temperature estimates back to 1300 AD. Figure 7 shows smoothed temperature anomalies with regard to the reference period 1989-2010.These anomalies are directly used to describe temperature variability for Nordenskiöldbreen since 1300 AD.
An accumulation history, derived from an ice core drilled at the top of the Lomonosovfonna plateau in 1997, provides estimates back to 1598 AD.An accumulation record back to 1715 AD, based on the similar ice core data set, is presented in Pohjola et al. (2002).The time series are extended to cover the period 1997-2010, using precipitation estimates from the Longyearbyen weather station (computing precipitation anomalies with regard to the period 1960-1996 and applying these to the Lomonosovfonna record since 1997).For the period 1300-1598, accumulation data were lacking and values are set to the mean accumulation over the period 1598-2010.The resulting smoothed record of accumulation anomalies since 1300 is shown in Fig. 7.
The ice flow model requires input of the surface mass balance and ice-surface temperature (at the snow-ice interface) as surface forcing.Both the mass balance and subsurface temperature depend on air temperature and accumulation.A coupled distributed energy balance-snow model, presented in Van Pelt et al. (2012), solves the surface energy budget and simulates melt water percolation, refreezing and runoff in the firn pack.The coupled model is applied to the grid of Nordenskiöldbreen to quantify the mass balance and subsurface temperature sensitivity by performing multiple model experiments with perturbed accumulation and air temperature over the period 1989-2010.Based on the previously discussed trends in the Longyearbyen instrumental record temperature perturbations are defined per season, e.g., a 1 K temperature perturbation in winter corresponds to a 0.47 K summer temperature perturbation.In this way, seasonality in climate change is taken into account.Seasonal differences in observed precipitation trends are much less pronounced (Førland et al., 2009) and therefore not accounted for in the coupled model experiments.Temporal variability in the climate forcing in the sensitivity runs is based on 3hourly climate input from a regional climate model RACMO (temperature and specific humidity) and the Longyearbyen weather station (precipitation and cloud cover).From distributed output fields of the sensitivity experiments altitudinal profiles of the annual mean mass balance and ice-surface temperature are computed for every climate perturbation.In order to enable selection of the mass balance and ice-surface temperature for every altitude and climate perturbation, curve fitting per height bin and further interpolation are performed to complete the diagrams, which are shown in Fig. 8.The diagrams in Fig. 8 serve as look-up tables for the mass balance and ice-surface temperature given the current altitude of the grid cell and the reconstructed temperature and precipitation anomaly.Hence, while running the ice dynamical model forward in time, mass balance and subsurface temperature values can be selected dynamically for every grid point by considering the current air temperature and precipitation anomaly (Fig. 7).With this approach, altitudinal feedbacks of the mass balance and subsurface temperature are implicitly taken into account.

Setup
Model iterations are performed on a grid with a 500 m resolution, containing a total of 3432 grid points.Every model www.the-cryosphere.net/7/987/2013/The Cryosphere, 7, 987-1006, 2013 iteration starts with an initialization run over the period 500-1300 AD, forced with constant-in-time, height-dependent climate input (Sect.4.2) to initialize the ice geometry and its thermodynamics (approaching a steady state in 1300).The prescribed surface climate during initialization corresponds to the mean climate for the period 1300-2007.Initialization is followed by a run over the period 1300-2007 with the model being forced with reconstructed time-and heightdependent climate input (Sect.4.2).The difference between modeled and observed surface heights in 2007, referred to as the surface misfit, is used to correct the bed after every iteration.
Outside the grid of Nordenskiöldbreen, the ice is allowed to evolve in a limited sense.A "force-to-thickness" mechanism is applied there, which involves prescription of an artificial mass balance required to force the surface height to approach present-day (DEM) values.The resulting slightly smoothed surface outside the grid resembles present-day surface heights and avoids severe time-stepping restrictions, i.e., very short model time steps, related to steep surface gradients.The grid of Nordenskiöldbreen is selected along the present-day ice divide, implying negligible ice flow across the grid boundary.The limited ice thickness evolution at the ice divide and its fixed position increase uncertainty in the reconstructed bed near the ice divide.Moreover, low surface slopes and small ice fluxes limit the bed to surface transfer amplitude and hence bed heights cannot meaningfully be determined.Therefore, Eq. ( 3) is not applied to adjust the bed at the estimated position of the ice divide and bed heights are determined by nearest-neighbor interpolation between iteratively adjusted bed heights on the grid and fixed bed heights outside the grid.Outside the Nordenskiöldbreen grid bed heights are determined by means of the perfect-plasticity assumption (Nye, 1951), assuming a constant applied stress equal to the mean determined from thickness and surface slope estimates at GPR bed data locations.A minimum slope of 0.025 is used to avoid unrealistically large ice thicknesses.Additionally, some measures are taken to prevent inconsistent corrections near the margins.This includes (1) interpolating the bed for grid points that are directly connected to ice-free grid cells, (2) interpolating the bed for modeled icefree grid cells that are ice covered in the reference surface from the DEM, (3) not adjusting the bed for grid cells that are ice free in the reference surface and ice covered at the end of a model iteration.Interpolation involves application of a spring metaphor to estimate unknown bed heights by interpolating between surrounding bed heights from inverse modeling (on the grid) and/or fixed bed heights (outside the grid).
In these experiments, ice is set to calve of at the presentday position of the calving front and hence the glacier front is not allowed to move beyond its present-day position.This implies some uncertainty in ice flow near the glacier A value for the relaxation factor K of 0.25 has been used in these experiments and leads to stable convergence of the inverse method.It is worth mentioning that climate feedbacks, like the mass balance-altitude feedback, enhance deviations of the surface height after a bedrock adjustment.This is one of the reasons to use a smaller K than in the synthetic experiments.
We start the first iteration with an initial bed generated through interpolation between known bed heights from the GPR data set and from the DEM (ice-free regions).Note that the data cover only a very small part of the grid, leading to a relatively smooth interpolated initial bed.For the tributary glaciers a constant initial ice thickness has been assumed to derive basal topography.The sensitivity of the reconstruction to the chosen initial bed is discussed in Sect.4.4.1.When validating the reconstructed bed, we discuss a possible bias introduced by including GPR data in the initial bed at the start of the inverse method.
As mentioned, modeling errors and surface data noise may lead to undesired corrections of the bed.A stopping criterion is needed to terminate the iterative procedure before the reconstructed bed gets too severely contaminated by unphysical variability.Synthetic experiments have indicated that with the number of iterations the rate of improve-ment of the surface misfit declines as smaller-scale bed features are recovered.The relative significance of modeling errors, as well as surface height data errors, hence increases with the number of iterations.Habermann et al. ( 2012) indicated that for problems involving both data and modeling errors, a so-called "L-curve criterion" can be used to define a reasonable point of termination (Aster et al., 2005).The L-curve criterion is based on the principle that the iterative method should stop when further reduction of the surface misfit comes at the expense of (unrealistically) large adjustments of the bed.In order to define the point of termination, the residual norm, i.e., the square-root of the summed squared surface misfit, is plotted versus the model norm in a log-log plot.The L2-model norm is used, which is the square-root of summed squared deviations of the reconstructed bed height relative to the initial bed.A bending point in this graph represents the critical point beyond which overfitting is likely to occur and hence defines a reasonable point to stop the iterative method (Fig. 10).Note that an alternative stopping criterion could be to stop the inverse procedure when the misfit between modeled and observed bed heights minimizes.With the L-curve criterion we illustrate a stopping principle that can be applied in presence and absence of (scarce) bed data.In other inverse modeling studies reconstructing basal conditions underneath glaciers, iterative procedures are often terminated when the surface misfit falls below a certain predefined tolerance ("discrepancy principle"; Maxwell et al., 2008;Arthern and Gudmundsson, 2010)  or when the rate of change of the residual norm drops below a predefined threshold ("recent-improvement principle"; Joughin et al., 2006;Sergienko et al., 2008) .In problems involving both modeling and data errors, it is generally hard to define the tolerance required in the discrepancy principle.On the other hand, the recent-improvement criterion depends strongly on the predefined value of the threshold, for which constraints are lacking.

Results
Next, we present results of application of the inverse method on Nordenskiöldbreen.In a first set of experiments (Sect.4.4.1), the approach is applied to a bed allowed to freely evolve everywhere, whereas in a second experiment (Sect.4.4.2),an alternative approach is applied in which the bed is fixed at locations where radar data are available.

Experiment I: unconstrained bed
First, we discuss results of the application of the inverse method to an unconstrained bed, i.e., the ice-covered bed is not constrained by any radar observations during the iterative method.The sensitivity of the reconstructed bed to modeling errors is studied by performing experiments with a perturbed (1) till strength, (2) initial bed and (3) surface climate forcing.Both spatial variability in the basal shear stress and bedrock topography exert an influence on the surface height and velocity distribution.The basal shear stress depends on the ice thickness, the presence of water, lubricating the bed, as well as the material till strength, which is known to be high for glaciers underlain by hard crystalline bedrock and much lower for soft-bedded glaciers underlain by sediment (Clarke, 2005).The till strength directly affects the amount of basal lubrication, which implies a strong influence on the ice thickness distribution.Since we lack observational constraints on the sediment distribution underneath Nordenskiöldbreen, we prescribe a spatially invariant material till strength in our experiments.Any surface height variability related to a non-uniform distribution of the material till strength is therefore not taken into account and is acknowledged as a source of uncertainty in our approach.
The sensitivity of the reconstructed bed to changes in the material till strength is investigated by performing runs with different values for the material till strength φ (Eq.2) with values ranging between 11 and 15 degrees.Snapshots of the reconstructed bed and the surface misfit for the φ = 13 • experiment after 1, 5, 10 and 23 iterations are shown in Fig. 9. Clearly, the surface misfit reduces with the number of iterations, indicating convergence of the approach.Starting from a relatively smooth initial bed, the reconstructed bed gradually acquires more smaller-scale features.In accordance with the GPR data (Fig. 6b) an overdeepening is apparent in the main flow with minimum bed heights around sea level.Fig. 10a shows the RMS surface misfit as a function of the number of iterations.A log-log plot of the smoothed residual norm (Fig. 10b) versus the model norm is used to define a point beyond which further iterations are likely to deteriorate the recovered bed.The bending point, corresponding to the point of maximum curvature is found here after 23 iterations.In all experiments discussed in this section the stopping point was found between 22 and 25 iterations.The surface misfit remains largest near the margins, where the modeled ice thickness is very sensitive to the applied surface forcing leading to uncertain bed corrections.This was the main reason to interpolate the bed in the immediate vicinity of the margin.
Figure 11 shows the reconstructed bed height and bed discrepancy for the sensitivity experiments with perturbed till strength (1st and 2nd column) and perturbed initial beds (3rd and 4th column).The bed discrepancy is computed relative to the final bed height in the φ = 13 • experiment.Perturbing the material till strength affects mean sliding velocities and the ice thickness.A lower till strength enhances basal velocities, thins the ice and hence raises the reconstructed bed in regions where sliding is significant.Largest bed discrepancies arise in locations where the surface slope is large, due to a higher sensitivity of sliding to changes in the basal shear stress.In cold-based higher regions the reconstructed bed is only slightly affected by a change in the till strength.
Regardless of the initial bed at the start of the first iteration, with the number of iterations the reconstructed bed always seems to converge to a similar bed profile, as illustrated in the synthetic experiments (Fig. 4).This is only approximately true in the current experiment, where a small misfit will always remain in the surface heights since the model is unable to produce an exact match of the surface heights in the DEM.The stopping criterion forces early termination of the inverse method, which implies that the recovered bed at the stopping point may still be influenced by the chosen initial bed profile.Figure 11 (column 3) illustrates that, when repeating the φ = 13 • experiment with an initial bed, which is uniformly lowered by 100 m, the remaining absolute bed discrepancy after 23 iterations is relatively small (on average 11 m).We performed an additional experiment to test the sensitivity of the reconstruction to an initial bed with perturbed spatial variability.For that purpose, an initial bed was constructed using the perfect assumption with a relatively high slope threshold (discussed in more detail later in this section).The initial bed is shown in Fig. 14c and lacks prominent features like the overdeepening along the main flow line.The reconstructed basal topography and the bed discrepancy are shown in Fig. 11 (column 4) and show similar variability as in the φ = 13 • experiment; for example, the overdeepening is apparent in the reconstruction.Remaining absolute deviations are on average 17 m (was 85 m at the start of the experiment).These results indicate convergence to a similar bed profile after many iterations for both experiments with perturbed initial beds.Nevertheless, due to preliminary stopping of the inverse procedure we can conclude that an initial bed estimate with a realistic spatial pattern and mean bed height does increase the accuracy of the final reconstructed bed.
We also study the sensitivity of the reconstructed bed to perturbations of the climate forcing.In a first sensitivity experiment, the φ = 13 • experiment is repeated with a constant-in-time surface forcing, using an air temperature and accumulation rate equivalent to the mean since 1300 AD.The resulting bed discrepancy with regard to the Table 1.Comparison of modeled and theoretically derived ice thicknesses with observations along the radar track (Fig. 6b).δH mean is the discrepancy between modeled and observed mean thickness (2nd column).The 3rd, 4th and 5th column present the root-mean-square error (RMSE), the mean absolute difference and the correlation coefficient (R), respectively.The mean observed ice thickness is 322 m. time-dependent forcing experiment is shown in Fig. 12 (1st column).This illustrates that applying a constant surface forcing raises the bed in higher regions (thinner ice) and lowers the bed in lower regions (thicker ice).The temperature time series (Fig. 7) reveal relatively high temperatures since 1900 AD, preceded by a relatively cold period since around 1600 AD.In the non-constant climate experiment the ice thickens significantly over the entire grid for the period 1600-1900 AD, due to a positive mass balance as well as a reduction of ice velocities, which is related to a drop in surface temperatures (advection of colder ice to the bed).Since around 1900 AD, relatively high temperatures induce more melting particularly in the lower regions, locally thinning the ice (Fig. 12).As a result of surface steepening, enhanced ice velocities gradually lead to thinning of higher regions as well.Since this dynamical thinning process is slow, relatively thick ice is still present in 2007 in the timedependent climate experiment.The relatively large and systematic bed discrepancy found when ignoring the time dependence of the surface forcing indicates the relevance of accounting for temporal variability when recovering basal topography or initializing ice flow models for forecasting experiments.In a second climate perturbation experiment, the air temperature forcing for the spin-up period 500-1300 AD has been lowered by 1 K.The bed misfit in Fig. 12 (2nd column) indicates a very small influence of the spin-up climate on the final reconstructed bed.Finally, we perturb the air temperature for the period 1300-2007 AD by −1 K (Fig. 12, 3rd column) and find a more pronounced effect on the bed misfit.
Next, the reconstructed beds in the sensitivity experiments are validated against observed bed heights, obtained using GPR along a track near the main flow line (Fig. 6b).A scatterplot of observed bed heights (a) and ice thicknesses (b) versus reconstructed values for the experiments with varying till strength is shown in Fig. 13.As mentioned, changing the material till strength affects the mean ice thickness and induces a shift in the mean bed height in areas where sliding is significant.Figure 13 illustrates that much of the variability in the observed bed heights and ice thickness is recovered by the model.Nevertheless, significant local errors remain, which can be ascribed to multiple factors.In fact, the bed height discrepancies in the scatterplot reflect the sum of all uncertainties involved in the inverse approach, including uncertainties in the model physics, climate forcing, surface data, bed data and the stopping criterion.Differences in mean thickness, the RMS thickness misfit, the mean thickness misfit and the correlation coefficient are computed for all sensitivity experiments and presented in Table 1.The lowest RMS misfit and a high correlation coefficient are found for the φ = 13 • experiment, which can hence be regarded as the best estimate.Perturbing the initial bed affects the correlation coefficient and RMS misfit only slightly, indicating robustness of the approach for different initial beds.Although best results are found for the φ = 13 • experiment, the perturbed initial bed experiments also indicate that a possible bias related to the inclusion of the GPR data in the initial bed topography at the start of the first iteration is small.Applying a constant surface forcing leads to a lower agreement between modeled and observed ice thicknesses, indicating the relevance of using a time-dependent climate forcing.Perturbing the spinup temperature has a very minimal effect on the quality of the reconstruction, whereas the effect of a perturbation of the 1300-2007 AD temperature on the reconstructed bed is more pronounced.It is noteworthy that the bed misfit in Fig. 13a changes sign at an altitude of ∼ 600 m a.s.l., which is likely to related to inaccuracy of the modeled position of the cold to temperate ice transition, which has a substantial impact on the ice thickness distribution.
Additionally, a comparison of observed ice thicknesses with theoretically derived values, based on the perfectplasticity assumption, is made.In the perfect-plasticity assumption (Nye, 1951), the local ice thickness H is directly computed from the surface slope s, while assuming a constant yield stress τ 0 : where g is the gravitational acceleration and ρ the density of ice.The perfect-plasticity method computes unrealistically large ice thicknesses when the surface slope approaches zero.Therefore, a minimum slope is defined, below which the ice thickness is not computed and determined by interpolation.Local yield stresses have been computed for all GPR grid points, using observed ice thickness and slope from the DEM.The mean of the local yield stresses is used as the constant yield stress in the perfect-plasticity method.Beds were reconstructed for three different slope thresholds (0.015, 0.02 and 0.04) and the resulting beds are shown in Fig. 14.A high slope threshold (Fig. 14c) implies that bed heights cannot be determined for a large part of the grid and are hence determined by interpolation.Conversely, with a low slope threshold (Fig. 14a) ice thicknesses are computed for a larger part of the grid, but the larger thicknesses are more uncertain.
Comparing the derived beds in Fig. 14 to the final bed in the φ = 13 • experiment (Fig. 9) reveals quite significant discrepancies.It should be noted that the perfect-plasticity assumption is a simplified form of the shallow ice approxima-tion (all shear is concentrated at the base) in which longitudinal stresses are ignored.The approach is therefore unlikely to perform well in areas of significant sliding.Correlating the perfectly-plastic beds against the radar observations (Table 1 Table 1. of the grid.Due to substantial refreezing in the accumulation zone ice-surface temperatures are substantially higher in the accumulation area than in the ablation zone (Fig. 8).
A transition from temperate to cold basal ice occurs around a bedrock height of 600 m a.s.l.In the ablation area, simulated subglacial temperatures depend mainly on cold ice advection, strain heating and frictional heating from sliding and cause the base to be temperate in areas with substantial ice fluxes and sliding, and cold in areas near the grounded margin.This is a common pattern for polythermal glaciers in Svalbard (Pettersson, 2004).

Experiment II: constrained bed
In the previous experiments, bed height data were not used to constrain the bed and were only used for validation purposes.
Next, we present an alternative approach in which the bed is held fixed in locations where bed heights are known from GPR observations.The iterative method then attempts to recover the bed for the remainder of the grid.This approach hence provides an alternative method that could be used when scarce bed data are available and the main goal is to estimate the bed for the grid not covered by the observations.We repeat experiment I with φ = 13 • , with the only difference being the fixed bed at radar grid points.The reconstructed bed, the remaining surface misfit and the difference between the reconstructed topography and the bed in experiment I are shown in Fig. 16.Comparing the surface misfit distribution in Fig. 16b to the pattern in Fig. 9 reveals that the misfit remains larger in regions close to the GPR data grid points (Fig. 6).This is expected as the bed is not freely evolving in these areas.The difference between the beds in experiments I and II (Fig. 16c) illustrates that bed heights of the grid points in the vicinity of the fixed bed grid points tend to compensate for the absence of adjustments of the fixed grid cells.For example, if a positive surface misfit occurs for a fixed grid point, the surrounding freely evolving grid points likely also experience a positive surface misfit and a negative bed correction.In the end, the cumulative negative bed correction will be larger than in experiment I, since the bed overcompensates for the lack of movement of the neighboring constrained grid cells.This explains the alternating bed discrepancy pattern, visible in Fig. 16c, with an average bed discrepancy close to zero.Far away from the fixed grid cells, the effect on the reconstructed bed is very limited.It is hard to say whether fixing the bed at bed data locations leads to an improved bed reconstruction in comparison with experiment I, since means to validate the unconstrained part of the reconstructed bed are lacking.One might argue that in experiment II the bed is more reliable at the GPR locations, but somewhat less reliable in the vicinity of the GPR locations, due to the aforementioned overcompensation of the bed.Both approaches provide useful means to initialize ice flow model for prognostic experiments, since the surface misfit at the end of the iterations is minimized.In absence of bed data, the approach in experiment I is the only way to reconstruct basal topography.

Discussion and conclusions
In this study, a simple bedrock recovery method is discussed and applied in synthetic experiments and in a real-world setting (Nordenskiöldbreen, Svalbard).The inverse method can effectively be used to reconstruct distributed beds when bed data is scarce or lacking.Requirements for the method include a distributed surface height dataset, an ice flow model  Idealized experiments with a synthetic ice geometry and a temporally constant surface forcing study the functionality and robustness of the inverse method in three dimensions.It is illustrated that local bed adjustments, despite their nonlocal surface expression, lead to effective reduction of the surface misfit and recovery of the reference bed with the number of iterations.The robustness of the method when changing the initial bed, relaxation factor, mass balance forcing and bump dimensions has been investigated.Main findings include that (1) results are shown to be insensitive to the initial bed, (2) an optimum value for the relaxation factor can be selected for which fastest stable convergence occurs and (3) larger bumps are more easily resolved than smaller ones.
In contrast to the synthetic experiments, application of the method to Nordenskiöldbreen involves dealing with timeand height-dependent climate input, leading to a non-steady final geometry, and dealing with modeling and data errors, implying a risk of overfitting.A time-and height-dependent surface forcing was constructed using reconstructed temperature and accumulation over the period 1300-2007 in combination with output from sensitivity experiments with a coupled distributed energy balance-snow model.An L-curve stopping criterion has been used to minimize contamination of the bed with undesired variability related to data and modeling errors.Bed height data along a radar track are used for validation of the modeled bed topography.Sensitivity experiments on an unconstrained bed all indicate convergence of the inverse iterative method.We find that (1) the chosen material till strength parameter significantly affects the mean modeled bed height in the sliding zone, (2) the initial bed has a small impact on the bed reconstruction, and (3) a constant-in-time surface forcing reduces the quality of the reconstructed bed.While all errors in the approach related to uncertainty in the model physics, termination of the correction method, surface height data and surface forcing, are absorbed in the reconstructed bed, we find correlations up to R = 0.89 when comparing reconstructed ice thicknesses with observations.The best estimate of the bed has a mean absolute deviation of 61 m (mean ice thickness 322 m).The remaining unexplained variability in bed heights can mainly be ascribed to uncertain model physics (and in particular the description of sliding).All together, the lack of spatial variability of the material till strength, inaccuracy of the modeled polythermal basal structure, uncertainty in the climate reconstruction, inaccuracy due to overfitting or underfitting in some areas and uncertainty in upscaling local radar height data onto a coarser grid contribute to errors in the bed.The approach in this study is shown to produce more accurate bed estimates for this glacier than obtained using the perfectplasticity assumption (using different minimum slope thresholds).In a final experiment, the bed was held fixed in locations where bed heights were known from the radar data.Here, the inverse method provides a useful tool to interpolate the bed between known bed heights.However, reconstructed bed heights in the immediate vicinity of fixed bed heights are likely overcompensating for the lack of movement of the fixed grid cells.
A better agreement between modeled and observed thicknesses could be obtained by more precise selection of a constant material till strength or including a spatial dependence of the till strength.Furthermore, extending the subglacial model with a water transport scheme would allow for detailed modeling of spatiotemporal variability in subglacial water pressures and basal lubrication.This is however beyond the main scope of this work, in which we aim to demonstrate convergence and applicability of the inverse method.The application to Nordenskiöldbreen sets an example of how the method can be used on a real glacier geometry.The functionality of the inverse approach does not depend on the accuracy of the surface forcing and the quality of the ice flow model, implying that the applicability of the approach is not constrained to a single model or glacier.Obviously, the accuracy of the reconstructed bed does depend on the choice of model and the reliability of the prescribed surface forcing.As this method can potentially be applied to estimate the thickness distribution of larger sets of glaciers and ice caps, e.g., to estimate ice volume, one could decide to use a computationally inexpensive ice flow model at the expense of detail in the reconstructed bed.Application to a set of ice masses could involve application of the inverse approach to individual glaciers, given a surface height data set, a climate forcing and an ice flow model, to estimate the summed volume contained in a set of individual ice masses.
A major advantage of the approach is that, in addition to reconstructing basal topography, it directly provides a means to spin-up ice flow models for forecasting experiments, since the method searches for a bed for which modeled surface heights approach DEM values and thermodynamics are initialized at the end of a time-dependent experiment.Alternatively, in case an accurate map of basal topography would already be available from observations, an approach as in Pollard and DeConto (2012) (iteratively using the surface misfit to update the sliding coefficient) can be used to constrain spatial variability in basal slipperiness.However, for many glaciers and ice caps, bedrock data are scarce or absent and the presented approach provides a reasonable method for spin-up.
The L-curve stopping criterion provides a reasonable means to terminate the iterative procedure before the recovered bed gets too severely polluted with undesirable noise.Whether adjustments of the bed are still an improvement depends on the relative significance of modeling errors in respect to the bed to surface transfer amplitude.Both the magnitude of modeling errors and the transfer amplitude may vary significantly from one place to another, implying that at the stopping point some areas are prone to overfitting, whereas in other places underfitting is likely to occur.In this study, a constant relaxation factor K has been used and no attempt has been made to apply weighted bed corrections, with stronger corrections in areas where bed corrections are more likely an improvement and weaker corrections in regions where overfitting is more likely to occur.Assigning spatially dependent values to K seems only feasible when a distribution of the magnitude of modeling errors as well as the bed to surface transfer amplitude can reasonably be defined, which is not straightforward.Examples of a possible spatial dependence of K could include gradually weaker corrections towards the margins, where modeling errors become more significant, and stronger corrections in high-sloping regions, where the bed to surface transfer amplitude is large.Since the inverse method performs best when the bed to surface transfer amplitude is large, the accuracy of the bed reconstruction is likely to be higher for glaciers with a larger slope and/or a higher slip ratio (in case a suitable flow model is available).
In future work the presented approach could be used as a spin-up method for prognostic experiments on Nordenskiöldbreen or other glaciers.Proper initialization avoids abrupt changes at the start of prognostic experiments related to the lack of spin-up of internal thermodynamics and the ice thickness distribution, which is inevitably needed to accurately simulate the near-future evolution of glaciers and ice caps.

Fig. 1 .
Fig. 1.The prescribed surface mass balance (a) and snow-ice surface temperature (b) forcing fields in the synthetic experiments.

Fig. 2 .
Fig. 2. Recovery of a single bump with an amplitude of 150 m in an otherwise flat bed.Panels (a) and (b) show reference surface and bed heights; (c) and (d) show snapshots of the surface height misfit and recovered bed height after n = 10, 20, 30 and 40 iterations.In (e) the gradual recovery of the bed is illustrated along the cross-section marked by the dashed line in (b); every line represents the bed after a certain iteration while the grayscale is a measure for the iteration number (low is light gray, high is black); the red line marks the reference bed height.Solid lines in (b), (c) and (d) mark a fixed position to illustrate movement of the recovered bump and surface misfit.

Fig. 3 .
Fig. 3. Recovery of a field of bumps with an amplitude of 150 m.The reference surface and bed are shown in (a) and (b).Panel (c) shows snapshots of the recovered bed height after n = 10, 20, 30 and 40 iterations.In (d) the gradual recovery of the bed is depicted with every line representing the bed after a certain iteration while the grayscale is a measure for the iteration number (low is light gray, high is black); the red line marks the reference bed height.

Fig. 4 .
Fig.4.The root-mean-square bed height misfit versus the number of iterations for experiments with varying relaxation factor K (green), initial bed (red) and balance gradient (blue).

Fig. 5 .
Fig. 5. Recovery of a bump field containing mixed bump sizes.The reference surface and bed are shown in (a) and (b).The reference bed contains a superposition of bumps with an amplitude of 150 m and a bump width varying by a factor 4. Panel (c) shows snapshots of the recovered bed height after n = 5, 20, 50 and 100 iterations.(d) shows the bed height misfit as a function of the number of iterations.

Fig. 6 .
Fig. 6.Observed surface heights from a DEM (a) and bed heights from radio-echo sounding (b).The black line marks the extent of Nordenskiöldbreen, connected along the ice divide to the Lomonosovfonna ice plateau.The location of mountain ridges Terrierfjellet, Ferrierfjellet and Minkinfjellet is indicated.The inset figure in (b) shows the Svalbard archipelago with the position of Svalbard Airport (SA) and Nordenskiöldbreen (NB).

Fig. 7 .
Fig. 7. Time series of reconstructed air temperature (red) and precipitation (blue) relative to the mean over the reference period 1989-2010.A 10 yr smoothing filter is applied.Precipitation anomalies are given in percentages, whereas temperature anomalies are absolute deviations.

Fig. 8 .
Fig. 8. Diagrams showing the surface mass balance sensitivity as a function of altitude for changes in air temperature (a) and precipitation (c).The sensitivity of the ice-surface temperature as a function of altitude for changes in air temperature and precipitation is shown in (b) and (d), respectively.Sensitivities were obtained from output of sensitivity experiments with a coupled energy balance-snow model, ran over the period 1989-2010.Temperature and precipitation anomalies are defined w.r.t the mean over the period 1989-2010.More details can be found in Van Pelt et al. (2012).

Fig. 10 .
Fig. 10.Plots showing the surface misfit vs. the number of iterations (a) and the residual norm vs. the model norm (b) for the φ = 13 • experiment.A bending point in the L-curve, used to define a stopping point for the iterative method, is marked in red.Five-point smoothing of the residual norm series facilitates selecting a bending point in (b).

Fig. 13 .
Fig. 13.Scatterplots showing the modeled bed height vs. the observed bed heights from radar (a) and modeled vs. measured ice thickness (b).Results are shown for sensitivity experiments with varying material till strength.The corresponding correlation coefficients, RMS misfit, bias and mean error are given inTable 1.

Fig. 15 .
Fig. 15.Distributed map of simulated ice velocities (a) and basal temperatures (b) in 2007 for the φ = 13 • experiment.At the locations of velocity observations, reported in Den Ouden et al. (2010), observed velocities (white) and simulated velocities (yellow) are indicated.Only basal temperatures higher than 270 K are shown.

Fig. 16 .
Fig. 16.Reconstructed bed with bed heights fixed when radar observations are available in panel (a).The corresponding remaining surface misfit is shown in panel (b).Panel (c) shows the bed discrepancy with regard to the basal topography in the unconstrained experiment with φ = 13 • .