Sensitivity of basal conditions in an inverse model: Vestfonna ice cap, Nordaustlandet/Svalbard

. The dynamics of Vestfonna ice cap (Svalbard) are dominated by fast-ﬂowing outlet glaciers. Its mass balance is poorly known and affected dynamically by these fast-ﬂowing outlet glaciers. Hence, it is a challenging target for ice ﬂow modeling. Precise knowledge of the basal conditions and implementation of a good sliding law are crucial for the modeling of this ice cap. Here we use the full-Stokes ﬁnite element code Elmer/Ice to model the 3-D ﬂow over the whole ice cap. We use a Robin inverse method to infer the basal friction from the surface velocities observed in 1995. Our results illustrate the importance of the basal friction parameter in re-producing observed velocity ﬁelds. We also show the importance of having variable basal friction as given by the inverse method to reproduce the velocity ﬁelds of each outlet glacier – a simple parametrization of basal friction cannot give re-alistic velocities in a forward model. We study the robustness and sensitivity of this method with respect to different parameters (mesh characteristics, ice temperature, errors in topographic and velocity data). The uncertainty in the observational parameters and input data proved to be sufﬁciently small as not to adversely affect the ﬁdelity of the model.


Introduction
Svalbard is an Arctic archipelago at around 80 • N. Nordaustlandet is the second largest island in this archipelago and has two major ice caps: Austfonna and Vestfonna. These ice caps represent one of the largest ice-covered areas in the Eurasian Arctic. The Vestfonna ice cap is the smaller of the two ice caps and covers about 2400 km 2 and ranges in altitude from sea level to 630 m a.s.l. It is composed of a small, but relatively smooth and flat area in the center (Pohjola et al., 2011b;Dowdeswell, 1986;Dowdeswell and Collin, 1990). A larger fraction is occupied by fast-flowing outlet glaciers (Figs. 1 and 2; Pohjola et al., 2011b). Some of these have been suggested to be of surge-type (Hagen et al., 1993). A large part of the total ice flux towards the margin is through these spatially limited areas of fast-flow units.
The role of glaciers as indicator of climate change, as well as their contribution to sea level rise, is widely acknowledged . The fast-flow units and surgetype glaciers are especially important in that context, but are however the most challenging glaciers to model. The high velocities are mainly achieved through basal motion (Clarke, 1987), as ice deformation in regions where the bed is at the Published by Copernicus Publications on behalf of the European Geosciences Union. (1:100 000,1990) and completed with the International Bathymetric Chart of the Arctic Ocean (IBCAO) (Jakobsson et al., 2008) in the sea. The Ahlmann Summit is located at the position A, and the ice core mentioned in the temperature discussion was drilled at the location C. The names of the different outlet glaciers are detailed in Fig. 2. The glacier outline is also based on the NPI maps. All data are used in the coordinate system WGS 1994, UTM zone 33 N. melting point limits velocities to a range of 1-10 m yr −1 . If unfrozen and wet basal conditions are maintained, stable fast-flow of outlet glaciers occurs. In the case of evolving basal conditions, acceleration and deceleration will occur in response to changes in the release of frictional heat and water production (Bougamont et al., 2011). The basal conditions are thus crucial for the observed flow regime as well as possible changes or even surges.
Here, we use the Full-Stokes model Elmer/Ice Zwinger and Moore, 2009) to study in detail the conditions at the bottom of Vestfonna (VSF) in 1995. Determining the optimal basal conditions from the glacier topography and the surface velocities is an inverse problem. The focus of this paper is the application of an inverse method and its technical aspects rather than an interpretation or analysis of the obtained distribution. A key motivation of this paper is that the sensitivity of the inverse method to errors in input data or choice of model parameters such as mesh resolution may play a critical role in the derived basal properties. Hence, we perform sensitivity tests.
Recently, several methods have been proposed to solve this particular inverse problem, for example Morlighem et al. (2010); Goldberg and Sergienko (2011); Gudmundsson and Raymond (2008); Raymond and Gudmundsson (2009); Pra- long and Gudmundsson (2011). Morlighem et al. (2010) study Pine Island Glacier, West Antarctica with a control method using an adjoint model. The spatial patterns of basal drag inferred are similar to those obtained with simpler models (Vieli and Payne, 2003;Joughin et al., 2009). Goldberg and Sergienko (2011) also use an adjoint model, while Gudmundsson and Raymond (2008); Raymond and Gudmundsson (2009); Pralong and Gudmundsson (2011) rely on a Bayesian method.
Here we use the variational approach with a Robin inverse method proposed by Arthern and Gudmundsson (2010); Gudmundsson (2011) and subsequently used in simulations of Variegated Glacier and the Greenland ice sheet (Jay-Allemand et al., 2011;Gillet-Chaulet et al., 2012). This algorithm can be implemented numerically without developing an adjoint model, but rather using the normal forward solution of the Stokes equations. We present the results obtained with the Robin method and compare them with results of a simpler approach: the use of mutually varying, but constant, sliding coefficients in different predefined areas. This emphasizes the role of the particular basal conditions at every outlet glacier of the ice cap.
In Sect. 2 the different data sets used in this study will be presented in detail. Section 3 describes the ice-flow model as well as the inverse method. In Sect. 4.1 we discuss some previous modeling efforts with simple sliding laws to motivate the necessity of using a more sophisticated sliding law and from Sects. 4.2 on the results obtained for the basal conditions in 1995 are presented and discussed.

Research site
Approximately 75 % of Nordaustlandet is covered by glaciers (Hagen et al., 1993), most of these of typical Svalbard polythermal type. Usually, marginal areas are cold and frozen to the bed, while thicker central parts are close to melting point. Varied topography of valleys, peaks, and resulting ice thickness changes produce the outlet glaciers characteristic of Vestfonna. The Vestfonna (VSF) ice cap was the focus of a recent IPY project (Pohjola et al., 2011b, and references therein) with several glaciological fieldwork programs focused on establishing the mass balance, identifying flow dynamics, remote sensing, snow chemistry and paleo ice sheet reconstruction. The Nordaustlandet ice caps were also the target of an expedition in the International Geophysical Year (IGY) during 1957-1958 when ice cap geometry, temperatures, mass balance and ice thickness were measured (Schytt, 1964). More recently, Austfonna has been visited by expeditions led by University of Oslo (Dunse et al., 2009(Dunse et al., , 2012(Dunse et al., , 2011. The research area is presented in Fig. 1. Here we focus on the smaller ice cap (VSF), which is more topographically variable and has relatively more valleys and outlet glaciers than the neighboring Austfonna (ASF) ice cap. A variety of fast outlets is concentrated in a small area on VSF (Fig. 2).

Topographic data
The rough surface topography ( Fig. 1) of the ice cap shows clear ridges in its central parts: the main ridge is oriented east-west and another ridge extends northwards from the eastern part of the main ridge. Between the ridges and the coast, there are pronounced fast-flowing outlet glaciers (Fig. 2). Most of the outlet glaciers are orientated northsouth. However, Franklinbreen, the largest outlet glacier, is situated on the northwestern side of the ice cap and some smaller ones drain into northeastern fjords. Some of them have been suggested to be surge-type glaciers (Hagen et al., 1993). For example, Braun et al. (2011) and Pohjola et al. (2011a) study the recent acceleration of the northwestern outlet glacier Franklinbreen. The bedrock data used are presented in Petterson et al. (2011). The data set consists of combined ground-based pulsed radar data collected in 2008-2009 and airborne radioecho soundings data acquired in 1983 and 1986 over the ice cap (Dowdeswell, 1986). The airborne radar did not give good bedrock depths in the thicker central parts of the ice cap due to scattering and absorption losses in the ice, while the ground-based radar was not deployed at the fast flowing outlet glaciers for safety reasons. The radar data were combined, and hence represent the ice cap at different periods; however, crossover analysis of the limited number of places where both surface and airborne radars registered bedrock in- Fig. 3. Digital elevation model for the bedrock elevation, overlaid are the ground-based (black) and airborne radar profiles (blue, gray: lack of bed return). In the land areas, this map is completed with the NPI map and in sea areas with the IBCAO data. However, these data were not used to establish the bedrock map from the individual profiles.
dicates that the difference in ice thickness change over the period 1986-2008 was estimated as no more than 7 m over the whole ice cap over the whole period, which is in agreement with a widespread thinning rate on VSF of 23 cm yr −1 during the period 1996-2002 (Bamber et al., 2005). The bedrock data are shown in Fig. 3. From this figure it is obvious that the southwestern part of the DEM is poorly supported by actual data (Petterson et al., 2011). There is also a poorly covered region at the northeast coast, but since the bed topography is smoother here, this introduces less error in the data interpolation. Also, the radar profiles on the southern coast are aligned with the direction of flow of the outlet glaciers, so that any bedrock ridges between the outlet glaciers may not be captured well in the DEM. The total areal fraction with a good spatial autocorrelation in the ice thickness interpolation was estimated to about 55 % (Petterson et al., 2011). The rest of the ice cap was mapped with lesser accuracy. The subglacial landscape undulates between elevations of −160 m and +410 m relative to sea level at an the average ice thickness of 186 m.

Surface velocities
Tandem phase ERS-1/2 1-day InSAR scenes acquired between December 1995 and January 1996 were used to calculate the ice surface velocity structure of Vestfonna. This data set is presented in Dowdeswell et al. (2008) and Pohjola et al. (2011a) in more detail; in the following it is referred to as "1995 velocities". Two independent look directions www.the-cryosphere.net/6/771/2012/ The Cryosphere, 6, 771-783, 2012 (ascending and descending orbits) have been used. Vertical components of the velocities have been neglected during the transformation to horizontal velocities. This assumption proved to be justified when comparing to land-based GPS data (Pohjola et al., 2011a). GPS data are available over large regions of the ice cap, excluding fast-flowing and hence too crevassed parts of outlet glaciers. In Fig. 2 the total horizontal surface velocity is shown. In addition, the directional information was also used in this study. In most cases errors are assumed to be smaller than about 2 cm per day (or 7 m per year). The displacement maps derived from satellite SAR data were geocoded to UTM projection at 100 m posting using the NPI DEM from 1990. We observe two very different flow regimes: a general pattern of very slow ice flow over the central area of the ice cap, with pronounced high velocity areas in the outlet valley glaciers.

Thermal regime
Svalbard has a maritime climate characterized by cooler summers and warmer winters than usually expected at this latitude. Mean monthly air temperatures on VSF during the short summer season do not exceed +3 • C, and in the winter averages are between −10, • C and −15 • C reaching minimum values between −25 • C and −40 • C depending on the year . At Kinnvika research station in the northeast of VSF winter monthly minimum temperatures were around −25 • C, and summer monthly maximum +4 • C (Pohjola et al., 2011b, and references therein). The mean annual surface temperature T surf is estimated using the mean annual air temperature at sea level T sea and an elevation lapse rate γ where S(x) is the surface elevation at the horizontal position x. We use a lapse rate of γ = 0.004 K m −1 (indicated by Schuler et al., 2007;Wadham et al., 2006;Wadham and Nuttal, 2002). T sea = −7.7 • C is estimated from two weather stations on Austfonna during 2004 to 2008 (Schuler, 2007) as well as four weather stations during 2008 and 2009 on Vestfonna .

Geothermal heat flux
Svalbard belongs to the Barents Sea shelf of the Eurasian Plate (Labrousse et al., 2008). The tectonic history of the archipelago is complex, and the rock outcrops of the coastal exposures of western Nordaustlandet formed during different geological eras. In southwestern Nordaustlandet Paleozoic strata rest on a Precambrian basement (Harland, 1997). Mesozoic dolerite intrusions dot the coastal exposures of the western part of the island. The basement of Svalbard consists of terraces assembled during the Caledonian orogeny in the Paleozoic era (Labrousse et al., 2008). The basal ice temperature depends on the geothermal heat flux, which varies with the bedrock type and age (Waddington, 1987). A heat flow value of 63 mW m −2 , representative of the post-Precambrian, non-orogenic tectonic region (Lee, 1970), is used for Vestfonna in this paper.

Temperature profile in the ice
In Svalbard glaciers are typically polythermal: at higher elevation on the ice cap, the release of latent heat when seasonal meltwater percolates and refreezes in the firn pack creates an additional heat source leading typically to temperate ice covered by cold firn. An example of such a typical accumulation area temperature profile under the influence of refreezing meltwater is shown in Fig. 4 (gray dots) (Watanabe et al., 2001). The ice core studied by Watanabe et al. (2001) was drilled at the east end of the summit ridge in June 1995 to a depth of 210 m where the over-all ice thickness is around 340 m (Petterson et al., 2011, Fig. 3). Theoretical temperature profiles based on the assumption of a steady state will lead to different shapes of the profile due to the differences in net mass balance in accumulation and ablation areas (Fig. 4, red and blue curve) (Paterson, 1994;Suter, 2002). In addition, in the accumulation area the profile is modified to become warmer in the upper part by the release of latent heat from refreezing of percolating melt water (Fig. 4, gray curve). At lower elevations in the ablation zone, the ice is usually below the freezing point, since solid ice prevents any percolation and refreezing of water. Thinner snow cover at lower elevations also allows efficient outgoing longwave radiative cooling of the glacier during the long The Cryosphere, 6, 771-783, 2012 www.the-cryosphere.net/6/771/2012/ Arctic winters. The profile for an equilibrium line is given in the same figure (Fig. 4, pink line). The temperature profile in VSF is poorly known, and in this paper we will use two extreme cases to investigate the sensitivity of the model to variations in the temperature profile. Naturally, the assumption of temperate ice throughout the whole ice body is an upper, even though unrealistic limit for the ice temperature. As simple estimation of the lower limit, we use the steady state equilibrium line temperature profile (pink line in Fig. 4); this scenario will be referred to as "cold scenario" in the rest of the paper. In the ablation area the real temperature profile will be slightly warmer as a consequence of advection of warmer ice from deeper ice layers. However, our approximation will be relatively close to a typical ablation area temperature profile -the area we are most interested in. In the accumulation area this "cold scenario" profile may be, strictly speaking, a colder limit only for the upper part of the ice column and slightly too warm in the lower part depending on the amount and penetration depth of latent heat. This temperature profile is applied in a de-coupled way from the velocity field.
By neglecting the advective heat transport, the equilibrium line temperature profile is determined by the surface temperature and the geothermal heat flux with φ the geothermal heat flux, κ the heat conductivity and D(x) the ice depth (distance to the surface at a given location x in the ice body). We use κ = 2.072 W K −1 m −1 (Ritz, 1987), φ = 63 mW m −2 , γ = 0.004 K m −1 and T sea = −7.7 • C (see above) as parameters determining T surf . Temperatures in Eq.
(2) are limited to values below pressure melting; T = 0 • C. We obtain bedrock temperatures between −8 • C near the coast and 0 • C in the thick central parts of the ice cap.
In summary, we assume that our idealized temperature profile is a lower limit for temperatures and will furthermore show that the results do not differ significantly from what we obtain with an isothermal ice body that corresponds to an (unrealistic) upper limit of possible temperatures.

Ice flow model
The modeling of glaciers and ice sheets has a large range of applications, such as ice core dating or simulating the role of large ice masses on the climate system. Whatever the objectives, numerical flow models seem to be the best way to capture the complexity of glacier evolution. The ice flow model solves the momentum and mass conservation equations for given initial and boundary conditions. The evolution of a glacier's surface geometry is then fully determined by the ice flow and the climatic forcing. For this it is crucial to understand the basal conditions and to use a suitable sliding law.  Paterson (1994) and other parameters used.

Forward ice flow model
A right-handed (x, y, z) Cartesian coordinate system is used. The surface and bedrock are defined as z = S(x, y, t) and z = B(x, y, t), respectively, but the bedrock will be assumed to be a fixed boundary, so that B(x, y, t) = B(x, y). Also, in this work only diagnostic simulations are performed, so that it is S(x, y, t) = S(x, y). The ice thickness is given by In what follows, the ice is considered as nonlinear, viscous, incompressible material. Due to ice incompressibility, the mass conservation equation reduces to Here, v x , v y , v z denote the components along the three directions of space of the velocity vector v. The stress tensor τ splits into a deviatoric part τ ′ and the isotropic pressure p where δ ij is the Kronecker symbol. The constitutive relation, which links deviatoric stresses to strain ratesǫ, is a power law, referred to as the Glen's flow law in glaciologẏ where A(T ) is the deformation rate factor (see Table 1 for adopted values), T the ice temperature (in this work variations of the pressure melting point are neglected and it is assumed to be at 0 • C) and τ ⋆ is the second invariant of the deviatoric stress tensor which is defined as The strain rate components,ǫ ij , are defined from the velocity components aṡ www.the-cryosphere.net/6/771/2012/ The Cryosphere, 6, 771-783, 2012 The force balance (quasi static equilibrium) in the three directions of space leads to the Stokes equations ∂σ x ∂x + ∂τ xy ∂y + ∂τ xz ∂z = 0 , where σ i stands for τ ii , g is the norm of the gravitational acceleration vector, g, and ρ is the glacier ice density (see Table 1).

Boundary conditions
At the upper boundary a stress-free boundary is assumed; the corresponding Neumann condition reads τ · n = 0, on S.
At the lower boundary we assume zero basal melting (v · n = 0) as well as a linear friction law (Weertman type sliding law, Robin type boundary condition (Greve and Blatter, 2009)). A Robin boundary condition is a linear combination of Dirichlet (in our case velocities) and Neumann (in our case stress) boundary conditions (Arthern and Gudmundsson, 2010).
where n and t are normal and tangential unit vectors, β is the basal friction parameter and v || and τ || are the basal velocity and stress components parallel to the bed. This basal friction parameter will be inferred in this study by an inverse method from measured surface velocities. We emphasize that β is the basal friction or drag parameter and not the commonly used sliding coefficient (i.e. higher values of β correspond to increased friction and less sliding).

Numerics
In the current study, this system of four partial differential equations with four independent unknowns (three components of the velocity vector and one for isotropic pressure) is solved numerically using the finite element method based on the code Elmer Zwinger and Moore, 2009 where scalar as well as vector-like test functions are introduced. The left-handed side term in the momentum equation has been integrated by parts and reformulated by applying Green's theorem. The deviatoric stress tensor τ ′ is expressed in terms of the strain rate tensorǫ by inversion of the power law (Eq. 5). The non-Newtonian stress strain relation introduces nonlinearities into the system implying the application of an iteration scheme. The system is solved applying the standard Galerkin method. Stabilization is obtained applying either the stabilized finite elements (Franca and Frey, 1992;Brezzi et al., 2004) or the residual-free bubble method (Baiocchi et al., 1993). Further discussion on the methods mentioned above is to be found in Gagliardini and Zwinger (2008).

Inverse model
Despite the available observations of surface properties, an important source of uncertainty in ice dynamics remains with the determination of basal drag coefficients. Here we use a method to determine this basal drag coefficient β from the more readily obtainable data set on the glacier surface (geometry and velocities) targeting the best fit between observed and modeled surface velocities.
The method is a generalization of the inverse method proposed by MacAyeal (1992MacAyeal ( , 1993 with the practical advantage that an already existing forward Stokes model is sufficient. It is a variational method, which relies on the minimization of a cost function measuring the mismatch between the modeled and measured surface velocities. Unlike earlier methods, the method used can be applied to a wider class of models like higher-order models or full-Stokes model similar to Maxwell et al. (2008). It can be applied to any spatially variable Robin coefficient.

The inverse method
We use the same approach as Arthern and Gudmundsson (2010); in the following we summarize its basic steps. The Stokes Eqs. (8) are solved iteratively with two different sets of boundary conditions. A Neumann condition at the surface represents the case where the surface velocity is freely determined from the ice geometry, basal friction parameter, etc., whereas the Dirichlet condition constrains the model according to the measured ice velocities. In the latter case, the upper surface stress field is allowed to vary in response. A cost function is defined and minimized when the Dirichlet condition on velocity and Neumann upper surface free-stress condition are both satisfied and hence determines when to stop the iteration. This Neumann-type problem is defined by Eqs. (3), (8) and the surface boundary condition (9). The associated Dirichlet-type problem is defined by the same equations, except that the Neumann upper surface condition (9) is replaced by the following Dirichlet condition (13) v hor (x) = v obs (x), ∀x ∈ S, The Cryosphere, 6, 771-783, 2012 www.the-cryosphere.net/6/771/2012/ where v hor (x) and v obs (x) stand respectively for the modeled and observed horizontal surface velocities. In z-direction, (τ · n) · e z = 0 is imposed all over the surface S, where e z is the unit vector along the vertical. These conditions are enforced for each location where surface velocity data are available. In absence of surface velocities, the Neumann Eq. (9) is applied for all three directions. The cost function, which expresses the mismatch between the two solutions for the velocity field with different boundary conditions on the upper surface S, is given by where the superscripts N and D refer to the solutions of the Neumann and Dirichlet problems, respectively. As shown by Arthern and Gudmundsson (2010), the inverse problem can be solved by minimizing this cost function. Furthermore, the Gâteaux derivative of this cost function with respect to β for a perturbation β ′ can be expressed as a The expression for the gradient in Eq. (15) is strictly valid only for linear rheologies, while the gradient for the nonlinear rheology of ice is unknown. However, Arthern and Gudmundsson (2010) showed that Eq. (15) also minimizes the cost function in the non-linear case of Glen's flow law, so we can use it here. Between iteration t and t + 1, β should be updated as β t+1 = β t − α β v D 2 − v N 2 , with a stepsize parameter α β . As we have large differences in the order of magnitude of the velocities between the interior and the margins, in this study an observational log-type law was used instead of this mathematical correct law. It proved to work well and leads to better convergence properties , with a step-size parameter ν.
We note that this log-type law is not purely empirical. Equalizing the original linear form of updating β and Eq. (16), it can be shown that this new rule corresponds to a particular (spatially varying) choice of α β , which is guaranteed to be positive, and vanishes as step-size parameter ν tends to zero. The minimization of the cost function J is obtained by the following iterative descent algorithm: 1. initial estimate of β, 2. solution of the Neumann problem, 3. solution of the Dirichlet problem, 4. update of β using an approximation for gradient of the cost function d β J , and 5. return to step 2 unless convergence is achieved.

Synthetic 3-D test case
Before applying the inverse method to real three-dimensional cases like the VSF ice cap or the Greenland ice sheet (Gillet-Chaulet et al., 2012), our implementation in the finite element code Elmer was thoroughly tested on a synthetic 3-D test case -the ISMIP experiment C (ice stream flow 1, Pattyn et al., 2008). Experiment C consists of an inclined ramp as bedrock topography, a constant ice thickness over the whole domain and a sinusoidal friction coefficient.
For different values of the domain length (sinusoidal frequency) and the slope of the bedrock, surface velocity distributions have been computed. Up to 2 % noise was added before using these velocity data as input as observed surface velocities in the inverse method. Different initial values for the basal drag coefficient were tested. The inverse method gave satisfactory and robust recovery of the initial sinusoidal friction distribution.

Technical aspects
The 1995 data set is sufficiently complete that small data gaps could be filled by interpolation, and the Dirichlet condition could be enforced everywhere. Iterations between the Neumann and the Dirichlet problem from an initial value β = 0.001 MPa yr m −1 over the whole domain are conducted. A β max threshold is introduced to prevent the model from producing unrealistically large values at underdetermined locations. That means that very small values and those close to β max do not have a physical meaning, but indicate that the model cannot adapt to fit the data. The value of β max was set to 1 MPa yr m −1 after testing that the obtained solution did not change compared to runs with 10 MPa yr m −1 or even higher values. The lower boundary is naturally zero.
The relative step power coefficient ν in the update of β (Eq. 16) was set to 3. In the simulations presented in this work, these iterations could be stopped after 15 iterations since the convergence of the cost function stagnated very quickly.

Finite Element Mesh
A critical factor in FEM is to discretize the 3-D glacier into a suitable mesh. Too large an element size results in missing crucial effects, while too small sizes result in unrealistically large computing resources (memory as well as computing time). We investigate in Sect. 4.3 effects of varying the horizontal and vertical resolutions. An ideal mesh should reflect the local physical variability scale. Our finite element mesh has been established using a fully automatic, adaptive, isotropic, surface remeshing procedure (Yams Frey, 2001). The mesh size is adjusted according to the 2nd derivations of the measured horizontal surface velocities. The obtained mesh (see fast-flow area and as coarse as 2.5 km in the slow-flow areas.
In the vertical direction we use 10 equidistant layers.

Application and sensitivity tests of the inversion of the basal drag force
To illustrate the importance of the inverse method solution, we firstly present a simpler approach. We then show the results for the basal drag coefficient using our best model setup and observational data followed by discussion of the sensitivity and robustness of the method as a function of various parameters: the mesh resolution, the surface geometry, the velocity data and temperature within the ice body.

Forward model using estimated sliding parameters
Studying the observed measured surface velocities, it is clear that sliding occurs on the VSF ice cap and features the fastflowing outlet glaciers. Attempts to model the ice cap with one and the same sliding law and parameter everywhere failed, and the model exhibited numerical instabilities. We do not show this elementary case here. Instead, we attempted to differentiate "by hand" areas with different coefficients in the sliding law. These two areas were characterized by a slow-flow region (observed velocities below 25 m yr −1 ) and a fast-flow region (velocities above this threshold). This leads to a stable model which reflects the existence of the outlet glaciers, but yet fails to reproduce the correct velocity pattern of those glaciers. In Fig. 6 (right) an example of this approach is shown. The poor velocity match was produced despite using sliding coefficients tuned to representative values in each outlet obtained from the inverse method (Sect. 4.2). The values are given in Table 2. There is ambiguity in arbitrarily deciding on the boundary between the two sliding areas and choosing the sliding coefficients. In contrast, the application of the inverse method does not require defining slow and fast-flow areas by hand or tuning of coefficients. Furthermore, this will lead to a correct flow pattern even inside the outlet glaciers.

Inverse model using 1995 data
The inverse method as described in Sect. 3 was applied to the VSF ice cap using the 1995 surface velocity data. The convergence was surprisingly fast, and after only 15 iterations the algorithm was stopped. The trend of the cost function (normalized to the value of the first iteration) throughout the minimization algorithm is presented in Fig. 7 on a logarithmic scale. The cost function does not decrease to zero since we are dealing with real data, and the absolute value has -if any -meaning only in terms of how well it is possible to fit the data with the model. The obtained basal drag coefficient β is presented in Fig. 8  (left). The blueish colors correspond to areas with little sliding, with the glacier is nearly frozen to the bed, while the reddish areas correspond to areas with very fast sliding. The overall pattern corresponds very well to that which may be expected from the surface topography and measured surface velocities. The limits of β are given (naturally) by 0 and (artificially) by 1 MPa yr m −1 . We examined (not shown) changes in β very close to zero (fast-flowing grid points). Introducing a lower threshold of β= 10 −6 MPa yr m −1 did not alter the resulting velocity field, showing that very low values of β do not affect the velocity field.
The velocities obtained with the model (see Fig. 6, left) are very close to the observed surface velocities. The histogram ( Fig. 9) quantifies the error between measured and modeled velocities and shows that, for the most grid points, the error is below 20 m per year. We recall that the error in the observation is about 7 m per year. Figure 10 provides, in addition to the cost function and the histogram, the spatial distribution of the relative error between the velocities ob-  (Table 2). the Dirichlet models. It is important to consider this spatial distribution because of the two distinct flow regimes. Over most of the glacier, the difference is less than 30 %. There is no noticeable difference in the quality of the results between fast and slow flow regions. Most of the grid points having higher errors are located in the fast-flowing outlets or are isolated points close to the margins, plausibly due to observational errors. The correlation plot between the measured and observed surface velocities is shown in Fig. 11 and illustrates once more the agreement over the whole range of velocities.
We observe some spurious values in some isolated points in the center of the ice cap, where we expect little or no sliding. There we obtain very low basal drag coefficients. These points are situated in the slow-flow areas right on the top of the ridges (Fig. 8, inset), where we expect the data to be erroneous since small errors in position may cause large errors in velocity vectors. Nevertheless, since the velocities are small here, the big error in β has little impact on the resulting velocity distribution and hence can be ignored.

Mesh resolution
A test run with 15 vertical layers led to identical results (not shown) to the one produced by 10 layers, meaning that in this case 10 are sufficient. A similar mesh as the previously used was established, but with coarser resolution in the fast-flow areas (minimum mesh size of 500 m compared with 250 m). The results are presented in Fig. 12. Looking at the example of Frazerbreen, it is clear that some features remain unresolved with a too coarse resolution. Also, the cost function is increased by 20 % with the coarser resolution showing that finer resolution allows for a better reproduction of the measured velocities. This illustrates the importance of the choice of an adaptive mesh and emphasizes the necessity of a fine resolution in the regions of interest. The changes in the results when switching from a homogeneous mesh to an adaptive mesh are even more striking (Fig. 12, right). The reason for this is probably not only the increase in resolution, but also the fact that the method is sensitive to local grid resolution.  9. Histogram of the difference between modeled and observed surface velocities. Note that the error on the measured surface velocities is 7 m yr −1 and that the last bin contains all entries with values higher than 280 m yr −1 .

Sensitivity of the method to the DEMs
To test the sensitivity of the method to errors in the surface DEM and also to verify if the difference in dates between surface DEM and velocity data epochs matter, we made test runs with slightly modified surface DEMs. The method seems to be rather stable with respect to such changes. Test runs for a surface DEM, which was modified by 15 times the current mean annual surface mass balance , were performed and no significant changes in the basal friction pattern could be observed. This modification of the DEM is much greater than what was observed on VSF between the 1950s (Schytt, 1964) and today. This modification was  bedrock DEM, since the obtained basal distribution still reflects strongly the surface velocity pattern.

Influence of errors in the velocity data
The error on the measured velocities is estimated to about 7 m per year by Pohjola et al. (2011a). A simulation with a random error between −10 m and +10 m per year on both components of the horizontal velocities was performed. The obtained distribution of the basal drag coefficient as well as the modeled surface velocities was not affected by this perturbation (not shown).

Influence of the ice temperature profile
The inverse method simulation was conducted for an isothermal temperate ice body. This represents an upper physical bound to ice deformation, though it is certainly unrealistic over much of the ice cap. Results presented in Fig. 8 (right) are rather close to those we obtained with our preferred temperature profile. The correlation between modeled surface velocities ( Fig. 13) with the cold and the warm temperature scenario confirms the good agreement of the results obtained with the different assumptions for the temperature. In addition in the same figure the correlation is also shown for the basal velocities. In fact, it is not surprising that the surface velocities are quite similar in both scenarios since the inverse method minimized the difference between observed and modeled surface velocities. But one could imagine that different assumptions for the temperature distribution lead to differences in the viscosity and consequently in differences in the distributions of basal velocities. But in our case, the total velocity is dominated by the sliding; the deformation plays hardly any role and we observe nearly the same correlation between velocity fields with different temperature profiles on the surface and at the basis. We can conclude that our earlier presented assumptions are good enough for the study of the basal conditions.

Discussion and conclusions
We have presented a 3-D application of the inverse method proposed by Arthern and Gudmundsson (2010). The distribution of the basal friction parameter underneath the VSF ice cap was inferred from the 1995 winter surface velocity distribution. We obtain a good representation of the observed velocities and can reproduce the fast-flowing outlet glaciers as well as the acceleration pattern within these glaciers. The average value of the error between measured and observed surface velocities is 18 m yr −1 compared to an error of 7 m yr −1 in the observed velocities. The improvement compared to the simpler models with constant sliding coefficients (as shown in Sect. 4.1) is clearly visible and shows the necessity of a good knowledge of the spatial distribution of the basal friction parameter β.
The robustness of the method with respect to spatial resolution and errors in the input data has been studied, and other technical aspects have been discussed. We have also shown that the obtained β distribution does not depend strongly on the temperature, justifying our approximation on the temperature profile in the ice body. The correlation between modeled surface velocities obtained with different assumption for the ice temperature is satisfactory. Too small values of β have little physical meaning, leading to the conclusion that, even though the spatial distribution of β is important, very small local variations in β are not visible in the velocity distribution. The knowledge of the basal conditions will be not only important in prognostic simulations, but also to improve our understanding of the mechanism of surging. Fig. 12. Comparison of the distribution of the basal drag coefficient with the usual (left), a coarser mesh (middle) and a regular mesh (right). The figure shows a zoom on two of the southern outlet glaciers. Fig. 13. Correlation between the modeled velocities for the two different temperature scenarios -for the surface velocities and respectively the basal velocities. Since the movement is dominated by the basal motion, there is hardly any difference between both velocity fields. All velocities are given in m yr −1 .