The Cryosphere A simple inverse method for the distribution of basal sliding coefficients under ice sheets , applied to Antarctica

Variations in intrinsic bed conditions that affect basal sliding, such as the distribution of deformable sediment versus hard bedrock, are important boundary conditions for large-scale ice-sheet models, but are hard to observe and remain largely uncertain below the modern Greenland and Antarctic ice sheets. Here a very simple model-based method is described for deducing the modern spatial distribution of basal sliding coefficients. The model is run forward in time, and the basal sliding coefficient at each grid point is periodically increased or decreased depending on whether the local ice surface elevation is too high or too low compared to observed in areas of unfrozen bed. The method considerably reduces large-scale errors in Antarctic ice elevation, from several 100s to several 10s of meters in most regions. Remaining ice elevation errors over mountain ranges such as the Transantarctics are further improved by parameterizing the possible effect of sub-grid topography in the basal sliding law, representing sliding in deep valleys. Results are compared with modern velocity data, and various sensitivity tests are described in Appendices.


Introduction
One major uncertainty in modeling continental ice sheets is the distribution of bed properties that determine the rate of sliding at the ice-bed interface.Basal sliding over hard beds is commonly described by a law relating sliding velocity u b to the basal shear stress τ b , such as where C is a basal sliding coefficient that depends on intrinsic bed properties such as small-scale roughness, and N is the effective pressure, i.e., ice overburden not supported by basal water pressure (Cuffey and Paterson, 2010).In models without an explicit hydrologic component, basal temperature is often used as a surrogate: where T b is the homologous basal temperature (relative to the pressure melt point), and f is zero for T b below some threshold, usually a few degrees to tenths of a degree C below freezing, ramping up to 1 at the melt point (e.g., Pattyn, 2010).In many large-scale models the same simple forms are also used for basal motion over deformable sediment, representing shearing within the sediment itself which is usually not modeled explicitly (cf.Howell and Siegert, 2000;Pollard andDeConto, 2003, 2007;Oerlemans and Nick, 2006).Values of C(x, y) used in large-scale models to represent hard rock vs. deformable sediment vary by many orders of magnitude, roughly 10 −10 to 10 −5 m a −1 Pa −2 for n = 2.In this paper we use essentially Eq. ( 2) with a weakly non-linear dependence on τ b (n = 2).Other types of sliding laws are discussed briefly below.The large-scale distributions of sediment vs. hard bed represented by C(x, y), and hydrology represented by T b or N , are probably the major sources of error in simulations of modern grounded Antarctic ice.Relevant data are sparse and/or indirect (e.g., Drewry, 1976;Studinger et al., 2001;Tulaczyk et al., 1998;Peters et al., 2006;Bingham and Siegert, 2009;Ferraccioli et al., 2009;Bell et al., 2011;Young et al., 2011).In comparison, other major factors are (i) better constrained by observations or experiments, such D. Pollard and R. M. DeConto: Distribution of basal sliding coefficients under ice sheets as surface and bed topography, surface velocity, mass balance, internal core temperature profiles and ice rheology, or (ii) arguably have lesser effects on ice geometry, such as non-uniform or anisotropic ice rheology, basal mass balance, geothermal heat flux, and neglect of longitudinal stresses in Shallow Ice Approximation (SIA) models, at least on large scales.Although Heberler et al. (2008) found climate variations and other parameters to be as important as sliding in modeling Fennoscandian ice-sheet geometry, their assumed range of sliding parameters was quite small compared to the potential sediment vs. hard-bed range.Our assertion is supported by Briggs et al.'s (2011) ensemble modeling of Antarctica, who found persistent errors that could not be reduced by adjusting model parameters with C(x, y) excluded.Whitehouse et al. (2012) found a similar strong sensitivity to basal sliding parameters in simulating the last Antarctic deglaciation.
Most previous paleoclimatic continental ice-sheet models have widespread (∼ 100s of m) errors in modern Antarctic surface elevations, with regional errors of 500 m or more (e.g., Ritz et al., 2001;Philippon et al., 2006;Pollard and De-Conto, 2009;Martin et al., 2011;Whitehouse et al., 2012).In this paper we take the view that: 1.Such large errors are likely to undermine paleoclimatic and future modeling applications, and it is important to reduce them considerably.For instance, simulations of past and future stability of the West Antarctic Ice Sheet (WAIS) could be seriously astray if modern ice thicknesses in its major drainage basins are in error by 500 m or more.
2. Much of these errors are due to erroneous prescription of intrinsic bed properties, C(x, y), and not so much to errors in the other basal-sliding terms, i.e., basal temperatures or hydrology (f (T b ) or N) in Eqs.(1) or (2).
In support of the latter point, we note that a number of largescale Antarctic models agree approximately on the locations of basal freezing vs. melting areas, given the same geothermal heat flux map (e.g., Pattyn, 2010, PD12).This paper describes a simple procedure to minimize errors in modern ice surface elevation by adjusting C(x, y).The resulting C(x, y) is mainly a model-derived estimate of the actual sediment vs. hard-rock distribution below Antarctica.The procedure tacitly assumes that all errors are due to unknown bed properties, and ignores any canceling errors due to imperfect basal temperatures and other model shortcomings.As ice models improve in the future and better observations of bed properties become available, canceling errors will hopefully be detectable, and model-derived C(x, y) maps can be validated.For now, we suggest that the risk of canceling errors is the lesser of two evils, worth taking in order to eliminate O(500 m) errors in modern surface elevation.
The inverse method does not depend on the exact form of the sliding law, and in principle is applicable to any smoothly varying relation between u b and τ b , perhaps even to multivalued and bounded-drag forms (Schoof, 2005;Gagliardini et al., 2007;Pimentel et al., 2010), as long as intrinsic bed quantities can be adjusted to increase or decrease u b for any given τ b and hydrologic conditions.This is not the case for plastic rheology, with τ b bounded by a given yield stress and no sliding for smaller τ b (e.g., Bueler and Brown, 2009), which is not amenable to this inverse method.
The method also does not depend on model details outside of the sliding law, but the model does need to be run long enough for the procedure to converge, on the order of at least 40 000 yr for continental Antarctica (see Sect. 6).In principle it would be preferable to use a high-resolution model with full-Stokes or higher-order dynamics to fully capture ice streaming regions and grounding-line zones, but that is currently infeasible for 40 000 yr time scales.Here, we use a coarse-grid model with a simpler hybrid treatment of longitudinal stresses, which allows long-term simulations while still producing reasonable streaming flow and grounding-line migration (Pollard and DeConto, 2007, 2009, 2012; henceforth PD07, PD09, PD12).

Previous inverse modeling, contrast with current method
A number of previous studies have attempted to deduce basal stresses or sliding coefficients under modern ice sheets and glaciers, pioneered by MacAyeal (1992MacAyeal ( , 1993) ) and adapted for instance by Vieli and Payne (2003) and Joughin et al. (2004).Those studies were applied to limited regions using the Shelfy Stream Approximation (SSA) equations appropriate for stretching flow, and relatively sophisticated control methods to rigorously account for the nonlocal nature of the dynamics.Recent approaches using similar variational or adjoint methods have been applied to the Pine Island/Thwaites Glacier areas (Joughin et al., 2009;Morlighem et al., 2010) and to continental Antarctica (M.Morlighem, personal communication, 2012).Also, Price et al. (2011) applied a simpler method to Greenland, and Le Brocq et al. (2009) linked a similar method with a basal hydrology model for West Antarctica.All these studies are based on fitting modeled ice velocities to observed surface or balance velocities, with ice thicknesses and elevations prescribed to modern observed.Other recent inverse methods (Arthern and Gudmunsson, 2010;Raymond Pralong and Gudmundsson, 2011;Jay-Allemand et al., 2011) are also based primarily on fitting modeled to observed ice velocities.
Here, a much cruder algorithm is used, fitting to surface elevations, not velocities.The ice-sheet model is run forward in time, and the basal sliding coefficient C is periodically adjusted at each point according to the local elevation error (ignoring the fact that the dynamical equations are nonlocal).This is continued iteratively until the model surface elevations converge to the best fit with those observed.It guarantees that the model will produce realistic modern ice thicknesses and elevations in subsequent runs with C(x, y) prescribed from the inversion procedure.The procedure makes the implicit assumption that modern Antarctic elevations and temperatures are close to equilibrium with modern climate, i.e., unequilibrated glacial isostatic adjustments remaining from the last deglaciation are small.Given that assumption and others, in principle the method should yield the same results as the other velocity-fitting methods above, because of the close relationship between balance velocities, ice thicknesses and surface mass balance (see Sects. 7, 8 and Appendix E).
Section 3 outlines the model formulation, and Sect. 4 describes the inversion procedure for C(x, y) and presents basic results.Section 5 improves results over mountain ranges using a modified sliding law that depends on sub-grid topographic relief.The rate of convergence of the inversion procedure is considered in Sect.6, and in Sect.7 model surface velocities are tested against a recent modern dataset (Rignot et al., 2011).The concluding Sect.8 qualitatively compares our C(x, y) results with earlier inverse studies, and discusses future directions.Appendices A to G present various topics including sensitivities to modern bedrock topography and other uncertain model inputs.

Model outline
The ice-sheet/shelf model used here is an updated version of that in PD07 and PD09.As described there, the ice dynamics are a heuristic combination of the scaled SIA and SSA equations for shearing and longitudinal stretching flow, respectively.A parameterization relating ice velocity across the grounding line to local ice thickness (Schoof, 2007) is imposed as an internal boundary-layer condition, so that grounding-line migration is simulated accurately without the need for very high (∼ 100 m) resolution (Schoof, 2007;Gladstone et al., 2010;Pattyn et al., 2012).A polar stereographic grid is used, with relatively coarse 40-km grid resolution that permits the numerous long runs needed for this paper; some tests at 20 km and 10 km are included in Appendix C and show that the results are essentially unchanged at the higher resolutions, including in ice stream areas.There are 10 unevenly spaced vertical layers, with standard treatments of ice temperature advection and vertical diffusion.The model is non-polythermal and has no explicit basal hydrology.Time steps range from 0.5 to 2 yr depending on resolution.All changes to the model since PD09 are described in PD12; changes that are particularly relevant for this paper are outlined here.
Surface mass balance is computed from observationally based datasets of modern climatological Antarctic precipita-tion and temperature (ALBMAP, Le Brocq et al., 2010; with accumulation from van de Berg et al., 2006).Simple lapserate corrections are made for elevation differences from modern, and a basic positive degree-day scheme is used for melt (PD12).The bedrock model is as in PD07 and PD09, with local asthenospheric relaxation towards isostacy, and nonlocal lithospheric elastic deformation; the ice-free equilibrium bed topography is derived from modern observed (Le Brocq et al., 2010), isostatically rebounded with all ice removed.A simple two-value pattern of geothermal heat flux is prescribed with 54.6 mW m −2 under East Antarctica and 70 mW m −2 under West Antarctica.This seems preferable to choosing one or another of available geothermal datasets (Shapiro and Ritzwoller, 2004;Fox Maule et al., 2005) which differ considerably from each other on regional scales; as noted in Appendix F, inverse results are insensitive to the choice of dataset.
Changes to the model physics since PD09 include a new parameterization of oceanic melting below floating ice, a calving scheme, and sub-grid fractional area at the edges of floating ice shelves.Other changes include a wider basal temperature ramp from no sliding to full sliding (−3 to 0 • C here, compared to −0.5 to 0 • C in PD09), and linear rather than log-linear weighting of the no-sliding and full-sliding coefficients.That is, in this paper whereas in PD09 In both cases, the weighting r is given by with the ramp-width temperature T r = −3 • C in this paper, and −0.5 • C in PD09.
Here, T b is the basal homologous temperature, C is the effective sliding coefficient used in the dynamics, C(x, y) is the sliding coefficient for T b = 0 • C, adjusted in the inversion procedure, and C froz = 10 −20 m a −1 Pa −2 (which is small enough to prevent any discernible sliding but is not exactly zero to avoid divide-by-zero exceptions in the numerics).It is unclear whether algebraic (Eq.3a) or geometric (Eq.3b) weighting of C froz vs. C(x, y), or something in between, is most realistic, and depends on how subgrid variations in basal stress are propagated upwards into the mean flow (Gudmunsson, 2003;Hindmarsh et al., 2006).Equation (3a) favors more sliding compared to Eq. (3b), and slightly improves results in the current model compared to PD09.
For many of the inverse runs in this paper, we are only concerned with grounded ice.Unless otherwise noted below, (i) grounding lines are constrained to modern observed locations, and (ii) ice fluxes across grounding lines and floating ice-shelf thicknesses are still predicted by the model, but a very crude "inversion" scheme is used for floating ice, whereby the sub-ice oceanic melt rate is adjusted locally at intervals so as to maintain floating thicknesses close to observed (similar to MacAyeal and Thomas, 1986;Joughin and Padman, 2003).This is not the focus of the paper, and is just an expedient to maintain realistic floating ice shelves while we concentrate on grounded ice.

Inversion procedure
The inversion procedure is very simple.The model (with grounding lines and floating ice constrained as described above) is run forward in time, starting from modern observed bed and ice surface elevations.As described in Sect.3, the model is forced by constant observed climatology (Le Brocq et al., 2010), with lapse-rate corrections for changing surface elevations.Basal sliding coefficients C(x, y) are initialized to the simple two-valued pattern shown in Fig. 1b (results do not depend on the initial C(x, y)).At intervals of t inv years, at each grid point with grounded ice, the local basal sliding coefficient C(x, y) in Eq. ( 3a) is adjusted by a multiplicative factor: where h s is the current ice surface elevation, h obs s is that observed, and h inv is a scaling constant.
During the inversion procedure, basal temperature is still allowed to influence sliding (Eqs.3a and 4).Adjustments to C(x, y) in Eq. ( 5) are only performed for non-frozen points with T b > −3 • C, and C(x, y) is constrained to not fall below 10 −10 m a −1 Pa −2 , representing hard rough bedrock.(If Eqs.3a and 4 are ignored and C(x, y) is allowed to fall to much smaller values mimicking freezing, then results become seriously degraded in subsequent non-inverse runs as shown in Appendix A).Adjusted C(x, y) values are also not allowed to exceed 10 −5 m a −1 Pa −2 , representing the slipperiest deformable sediment.
These minimum and maximum limits of 10 −10 and 10 −5 for C(x, y) are themselves quite tightly constrained by model behavior.With a higher minimum, e.g., 10 −9 , basal flow is too fast in many regions and the modeled East Antarctic Ice Sheet becomes generally too thin; much lower values of ∼ 10 −12 are presumably unrealistic because they imply almost no basal sliding even for the largest τ b .With a lower maximum, e.g., 10 −6 , ice elevations in the Siple Coast region become too high with insufficient sliding in streams; if 10 −4 is used, the model is numerically unstable.
In Eq. ( 5), the interval between adjustments t inv , the scaling constant h inv , and the adjustment-factor limits (10 −1.5 to 10 1.5 , i.e., ∼ .03 to 30) are chosen to avoid overshoots and produce reasonably efficient convergence of the inver- sion procedure.For most results in this paper, we used t inv = 5000 yr and h inv = 500 m, and the runs all converged to essentially invariant states after ∼ 200 000 to 400 000 years (see animations in Supplement).We later experimented with other choices of t inv and h inv , and found that spinup times can be reduced considerably to ∼ 40 000 yr (see Sect. 6).

Basic results
Basic results of the inversion procedure applied to continental Antarctica are shown in Fig. 1.For purposes of comparison, the first column of panels in Fig. 1a-c shows noninverse results with a very simple two-value prescription of basal coefficients, i.e., a hard-rock value where ice-free isostatically rebounded modern bed elevations are above sea level, and a deformable-sediment value where they are below (Studinger et al., 2001;PD09;cf. Whitehouse et al., 2012).
Here, these values are 3 × 10 −9 and 3 × 10 −8 m a −1 Pa −2 , respectively, which produce more or less the smallest overall surface elevation errors, but much the same results are obtained with other pairs such as 10 −10 and 10 −5 .As shown in Fig. 1a, the departures from modern ice surface elevations are large, 500 m to 1000 m in places, and are typical of those in previous large-scale long-term Antarctic modeling mentioned above.
Results with the inversion procedure are shown in the second column of Fig. 1d-f.Surface elevation errors are generally much smaller compared to Fig. 1a.The improvement is shown quantitatively by histograms of error magnitudes in Fig. 2, which cluster around a mode of ∼ 20 m with the inversion procedure compared to ∼ 300 m with the two-valued C. Similarly, Table 1 shows mean elevation error magnitudes for most experiments in this paper.For the various inversion runs that are bona fide modern simulations, they range from 42 to 72 m, compared to 235 m for the two-valued C.

Inversion results with sub-grid topographic influence
Although surface elevations produced by the inversion are generally within a few 10s of m of observed, they are still too high by several hundred meters over most of the Transantarctics and some other mountain ranges (Fig. 1d).There the model's basal temperatures are uniformly frozen, and the inversion procedure cannot compensate for the hindrance to cross-range flow (cf.Kerr and Huybrechts, 1999).To some extent, frozen basal temperatures are expected over mountain ranges because the thinner ice provides relatively little insulation from cold surface temperatures (Pattyn, 2010).However, there may still be significant basal sliding in deep and warmer valleys not resolved by the coarse grids used here.We attempt to parameterize this sub-grid process by modifying the width of the basal-temperature ramp T r in Eq. ( 4) as a function of sub-grid topographic variations.The constant  1a.Green: using inverse method as in Fig. 1d.Blue: using inverse method with sub-grid topographic influence as in Fig. 3d.
value T r = −3 • C used above is replaced by where SA is the mean sub-grid slope amplitude computed by averaging the bed slopes in the 5-km ALBMAP dataset (Le Brocq et al., 2010) within each model grid box.h eq b is the ice-free isostatically rebounded (and 9-point smoothed) bed elevation on the coarse model grid, discussed below.SA was also used by Marshall et al. (1996) in another context.Whitehouse et al. (2012) apply a similar increase in sliding coefficient over mountainous terrain, for much the same reasons.
SA is typically ∼ .02or less in plains, and ∼ .03 to .05 or more in mountain ranges where T r can typically be ∼ −15 • C or colder.However, around the Gamburtsevs and also the extreme southern Transantarctics (∼ 120 to 180 • W), SA seems to be anomalously low, presumably due to the sparcity of BEDMAP data lines (Lythe et al., 2001).Hence, the second term in Eq. ( 6) uses the grid-scale elevation h eq b as a surrogate.h eq b is less than 1700 m nearly everywhere except over those two regions, where it is typically ∼ 2000 m or more; so again, T r can typically be ∼ −15 • C or colder.The use of h eq b is based only on the assumption that very high regions also have high sub-grid variability; future improvements in topographic coverage will probably allow just the first term in Eq. ( 6) to be used (Bo et al., 2009;Bell et al., 2011;Young et al., 2011).
Using Eq. ( 6), results with the inversion procedure are improved (Fig. 3d-f; second column).Now, more basal sliding over mountain ranges occurs despite frozen basal temperatures, reducing too-high surface errors.Remaining elevation errors in Fig. 3d are mostly < 50 m, except in a few small patches over mountains.As noted above, an important feature of the inversion procedure is that if the resulting C(x, y) distribution is prescribed in a subsequent non-inverse run with nothing else changed in the model, the results including surface elevations remain unchanged; for instance, with the modified sliding dependence in Eq. ( 6) and prescribed C(x, y) from Fig. 3e, the results would be exactly as in Fig. 3d and f.
In all simulations to this point, the focus has been on replicating modern Antarctic grounded ice, so grounding lines have been fixed at their present positions, and floating ice thicknesses have been artificially constrained to remain close to observed.Fig. 3g-i (third column) shows that when these constraints are relaxed and the complete icesheet-shelf model is integrated forward to equilibrium with C(x, y) prescribed from Fig. 3e (i.e., a non-inverse run with free grounding lines and shelves, and all other physics and forcing as above), there is only a little degradation of ice elevations, and modern errors remain less than ∼ 50 m nearly everywhere.The Ronne grounding line recedes slightly too much in the interior Weddell embayment, causing negative ice elevation errors there (Fig. 3g); also, the model fails to simulate the George VI Sound and ice shelf between Alexander Island and the western Peninsula, causing positive el-evation errors.However, overall the unconstrained model, with prescribed basal sliding coefficients C(x, y) from the inverse method with sub-grid topographic influence, produces a modern ice distribution very close to observed.

Convergence rate
All of the inverse runs described above use t inv = 5000 yr (interval between adjustments) and h inv = 500 m (scaling constant) in Eq. ( 5), and converge to a nearly invariant state after about 200 000 to 400 000 yr.Some slow, small regionalscale variations (few 10s of m in ice elevation, few 10 000s of yr timescale) continue indefinitely, but with very small effects on the deduced C(x, y) patterns.Two animations of a typical inversion run are provided as Supplement, showing model-minus-observed surface elevations and log 10 (C(x, y)) every 5000 yr through the 400 000-yr integration.
Spin-up times of ∼ 200 000 yr are feasible with our current hybrid model, but would be a serious impediment for more CPU-intensive higher-order and full-Stokes models.We have experimented with other choices of t inv and h inv to see if spin-up times can be reduced.The result shown in Fig. 4 are encouraging; for instance, with t inv = 500 and h inv = 1000, convergence of the mean elevation error takes on the order of 40 000 yr, a ∼ 5-fold speedup.All the runs shown in Fig. 4 asymptote to very nearly the same C(x, y) distribution by the time that the mean elevation error drops to ∼ 50 m.
Reasonable choices of t inv range from a few hundred to a few thousand years, which allow ice elevations to change appreciably between iterations but are not too large to cause overshoots.Reasonable choices of h inv are ∼ 50 m to 500 m, producing substantial changes in C(x, y) for elevation-error sizes of interest.We have not found any combinations of t inv and h inv that produce shorter spin-up times than 40 000 yr, or much different curves than those in Fig. 4.

Comparison with observed velocities
As discussed in Sect.2, in contrast to previous inverse modeling, the current method is based on fitting to observed ice elevations, not velocities.Nevertheless, our surface velocities should agree with those observed, given a number of conditions: 1. the modern surface mass balance dataset used here is realistic, A new all-Antarctic dataset of surface velocities (Rignot et al., 2011) provides the opportunity to test this, as shown in Fig. 5 where the dataset (900-m spacing) has been regridded by simple area-averaging to the model's 20-km grid (see Appendix C for other results at this resolution).Quantitative comparison is hindered by the fine scale and sharp gradients of many features in the dataset such as numerous outlet glaciers around the coast, many of which are barely resolved by the model and may be slightly displaced to one side or the other.Model speeds in the flanks around most coastlines are generally too fast, both in outlet glaciers and in the slower flows between them.The model's marginal ice thicknesses are generally close to observed (Fig. C1d), so the discrepancy might be caused by too much snowfall near the coasts, or too much internal deformation compared to sliding.The biggest single velocity error in Fig. 5c is due to the Kamb Ice Stream (Ice Stream C) on the Siple Coast, which stagnated about 150 yr ago (Hulbe and Fahnestock, 2007) but in the model is flowing at velocities comparable to the other active Ross ice streams.

Summary and discussion
A simple inverse method of adjusting basal sliding coefficients to minimize modern ice surface elevation errors (Eq.5) drastically reduces these errors in our continentalscale Antarctic model.Unlike the more sophisticated control or adjoint methods used in previous inverse studies, the method is local and ignores the spatial connectivity of ice dynamics.With basal temperatures included in the sliding parameterization during the inversion procedure, realistic ice elevations are maintained in subsequent non-inverse simulations with prescribed sliding coefficients.
To further reduce small patches of surface elevation errors (∼ 100s of m) over mountain ranges such as the Transantarctics, the influence of basal temperature on sliding had to be modified to include sub-grid topographic variations (Eq.6).Without this, frozen basal temperatures hinder sliding too much across major mountain ranges, causing local ice to be too thick.This parameterization of basal flow in deep sub-grid valleys may be tested in future work by higherresolution models (cf.Egholm et al., 2011).
Clearly there is a danger of canceling errors, i.e., the deduced sliding coefficients C(x, y) may not be real, but instead might be compensating for errors in the model physics or other input datasets.Our C(x, y) probably do represent a combination of (i) model errors, (ii) omitted physical variables affecting sliding, such as overburden pressure or hydrologic regime, and (iii) all intrinsic bed properties that affect sliding, not just sediment vs. hard rock but also small-scale roughness.In the future it may become possible to sort this out with better (higher-order, higher-resolution) models and new improved datasets, and achieve convergence between models and data on real Antarctic bed conditions.For now, we suggest that the danger of canceling errors is preferable to www.the-cryosphere.net/6/953/2012/The Cryosphere, 6, 953-971, 2012 using ice models with very large (>∼ 500 m) elevation biases to study important problems such as past and future stability of WAIS.Some first steps in checking for canceling errors are taken in appendices.Appendix B establishes rough bounds on the internal-flow enhancement factor E, outside of which the basal inversion cannot maintain a realistic ice distribution.Appendix F tests the sensitivity to uncertainties in some prescribed input fields -geothermal heat flux, surface accumulation, and unloaded bedrock -and finds that the deduced distribution of basal sliding coefficients C(x, y) is robust on large scales.
In principle our results can be compared with previous inverse studies using relatively sophisticated control or adjoint methods mentioned in Sect. 2. Although these studies fitted to observed velocities as opposed to surface elevations, the resulting basal maps should be similar due to the interdependence of equilibrium velocities, ice thicknesses and surface mass balance.Published Antarctic studies have targeted limited regions at higher resolution, rather than continental Antarctica as in this study.Also, most have deduced basal stresses or used different sliding laws than ours, so comparisons with our basal coefficient maps can only be qualitative.Joughin et al. (2009;their Figs. 5a, 6a) and Morlighem et al. (2010; their Fig.2a-c) deduced basal stresses in the Thwaites (TG) and Pine Island Glacier (PIG) areas.Comparing their maps with our highest-resolution C(x, y) pattern in Fig. C1h, the overall agreement is sporadic at best, although there are some zeroth-order points of agreement: (i) upstream of the PIG grounding line, there is a few 10s of km strip with higher basal stress and lower C values, then a ∼ 100-km-long broadening zone with low stresses and high C values; (ii) ∼ 100 to 200 km upstream (southward) from the TG grounding line, there are transverse strips of alternating high/low stresses and high/low C values, similar to more numerous strips in Joughin et al. (2009).
Directions for future work include combining the basal inverse method with statistical ensemble techniques involving other model parameters (Hebeler et al., 2008;Stone et al., 2010;Applegate et al., 2012;Briggs et al. 2011;Tarasov et al. 2012); this could test our assertion that basal properties are one of the largest sources of uncertainty in ice-sheet modeling.The inverse method could also be combined with basal hydrologic modeling following Le Brocq et al. (2009).The assumption of modern quasi-equilibrium, especially for internal and basal temperatures, can be tested by integrating through the last few 10 000s of yr (although we have previously found that modern ice thicknesses and groundingline locations are much the same in modern equilibrated vs. transient runs, not shown).Whitehouse et al. (2012) have recently developed a deglacial model of Antarctica from 20 ka to the present, and performed an extensive exploration of model parameter space, constraining or comparing to much of the available data on ice extents, thicknesses and elevations over this period.Analogously to here, they manually adjusted basal sliding coefficients in some regions.Perhaps the automated inversion procedure could be combined with transient deglacial simulations, which would allow basal sliding coefficients to be adjusted on the continental shelves, and could include other constraints such as relative sea level records (Bassett et al., 2007;Briggs et al., 2011), modern uplift rates (Ivins and James, 2005;Thomas et al., 2011), and adjustments to equilibrium bed topography (Gomez et al., 2010;Raymond Pralong and Gudmundsson, 2011).
This paper only addresses the modern distribution of Antarctic sediment vs. bedrock under current grounded ice, which has largely been shaped by erosion, transport and deposition by glaciers and ice sheets over the last ∼ 34 million yr (Jamieson and Sugden, 2008).If sediment distribution and bed properties can realistically be deduced in modern ice-sheet simulations, that will help in developing coupled ice-sheet-sediment-bedrock models (Jamieson et al., 2010) aimed at the long-term landscape evolution of Antarctica.

Appendix A Inversion without basal temperature influence
In Sect.4.1 we emphasized that the inversion procedure needs to account for the effect of basal temperature on sliding, i.e., Eqs. ( 3a) and ( 4) should be applied during the inverse run, as they are in non-inverse model simulations.If this is not done, very small elevation errors can be achieved in the inversion procedure, but results become seriously degraded in subsequent non-inverse runs.This caveat is illustrated here.
Figure A1a-c shows results with Eqs. ( 3a) and ( 4) turned off during the inverse procedure, i.e., with r set to 1 in Eq. (3a) (and without sub-grid topographic influence).Adjusted C(x, y) values are allowed to fall far below 10 −10 m a −1 Pa −2 to 10 −20 , which renders basal sliding insignificant even for the highest shear stresses so that the inversion procedure effectively predicts its own "frozen" areas.Surface elevation errors fall to very small values during the inversion procedure (< 50 m almost everywhere in Fig. A1a, mean of 23 m in Table 1), even without sub-grid topographic influence.
The deduced C(x, y) map (Fig. A1b) has deformablesediment values (10 −5 ) in streaming ice regions such as the Siple coast and other marginal outlet channels, as expected.C(x, y) values less than about 10 −12 , shown as purple in Fig. A1b, allow essentially no basal motion and imply that the bed should be frozen there.Note that the Transantarctics and most of the Dronning Maud Land mountains do not have these low values, and have high sediment-like values in some channels, i.e., the inversion procedure requires that ice slides easily over these mountain ranges (cf.Sect.5).
The purple regions in Fig. A1b can be compared with the frozen-bed areas in Fig. A1c (T b < −3 • C) predicted by the model thermodynamics.Although similar on the broadest scales, there are substantial regional differences, in particular over the Transantarctics which have quite high sliding coefficients in Fig. A1b but frozen basal temperatures in Fig. A1c.One can anticipate that surface elevation errors will worsen in these regions if the model is run in non-inverse mode with C(x, y) prescribed from Fig. A1b, and with T b allowed to influence sliding in Eq. ( 4).
The second column of Fig. A1d-f shows just that.The deduced basal coefficient map in Fig. A1b is first filled by replacing patches of C < 10 −10 with nearest-neighbor values ≥ 10 −10 .The resulting C(x, y) map (Fig. A1e) is prescribed in a non-inverse run in which basal sliding is influenced by basal temperature T b (Eqs.3a and 4).As expected, errors in ice surface elevations in Fig. A1d are much worse than in Fig. A1a.For instance, over the Transantarctics the model predicts a frozen base and no sliding, and ice surface elevations are ∼ 500 m or more too high.
To improve matters, one strategy would be to attempt to tune the model's thermodynamics to replicate the implied frozen regions of the inversion procedure (purple patches in Fig. A1b), perhaps by including new hydrologic processes (Bell et al., 2011).Although that might turn out to be a useful future direction, we note that the basal temperature (T b ) patterns cannot be made consistent with Fig. A1b by "easy" changes such as alternate geothermal heat flux maps (PD12) and different parameterizations of ice conductivity and specific heat (not shown).Furthermore, the modeled T b is much the same in all cases here (bottom panels of most figures), and mostly agrees on large scales with several other models (e.g., Pattyn, 2010) -frozen on buried mountain ranges and slow-flowing marginal flanks where thinner ice provides little insulation from the cold surface, and melting elsewhere, especially in fast-flowing streams and focused outlet channels.Hence, we have not attempted to tune T b , and have adopted the alternate inverse strategy described in the main paper with basal temperatures affecting sliding.

Appendix B Bounding of enhancement coefficient
As noted above, one concern with any inverse procedure for basal conditions is the possibility of canceling errors in other parts of the model or input fields.One example is potential errors in the internal deformation flow.Here we attempt to show this is not the case, at least to zeroth order, by performing a tuning exercise to crudely constrain the internal-flow enhancement factor E. (This factor multiplies all strain rates ∂u/∂z in Glen's law within the SIA dynamics.Enhancement factors are commonly used to improve large-scale geometries in ice-sheet models; e.g., Ritz et al., 2001).
Each column of panels in Fig. B1 shows results for a given E, from E = 0.1 to E = 8.With smaller E (0.1 to 0.5, first two columns), there is very little internal flow, and the inversion procedure attempts to compensate by assigning higher values of C(x, y) everywhere -but it cannot where the bed is frozen, so ice elevations in some regions become too high.The third column with E = 1 gives the best overall results (same as in Fig. 1d-f).For larger values of E (2 to 8, last three columns), there is too much internal flow, often exceeding the balance velocities implied by the surface mass balance.The inversion procedure cannot completely compensate even by setting C(x, y) to the smallest (hard rock, 10 −10 ) value, and surface elevations are generally much too low.We conclude that E = 1 is the most realistic internalflow enhancement factor for this model, and have used that value for all other simulations in this paper.Rignot et al. (2011) performed a similar estimate of the deformational factor for modern Antarctica by fitting to observed surface velocities in divide regions where internal deformation is expected to dominate ice motion.As they discuss, most internal deformation occurs near the base in ice within a few • C of the melt point.Their deduced creep parameter value of 9 × 10 −25 s −1 Pa −3 agrees well with ours, which for E = 1 and homologous temperature of −5 • C (mentioned in their Supporting Material) is 14.6 × 10 −25 s −1 Pa −3 (PD12).
One could conceivably attempt to further refine the model's internal flow physics, perhaps with spatially varying E or anisotropic rheology (e.g., Wang and Warner, 1999;Graversen et al., 2011), but these refinements would likely be beyond the scope of the present approach.As discussed in the introduction, one assumption of this paper is that largescale ice elevation errors caused by deficiencies in internal flow are minor compared to those caused by unrealistic bed sliding.

Model resolution
To better compare with earlier studies, and also to test for grid-size dependence within our model, Fig. C1 shows results at somewhat higher resolutions of 20 and 10 km.Compared with the corresponding 40-km results in Fig. 3d-f, some finer details emerge, but the large-scale pattern of basal coefficients and the amplitude of surface elevation errors remain essentially the same at the higher resolutions.This is true even for ice streaming regions such as the Siple coast, where 40-km resolution is not expected to resolve individual ice streams.The proto-streaming in the 40-km grid requires the same sliding coefficients as in the finer grids to produce the same regional ice thicknesses.This insensitivity to model resolution is reassuring, but may simply be due to intrinsic limitations of inverse methods to resolve smaller-scale basal sliding features (Gudmundsson and Raymond, 2008).

Appendix D Model bed elevation errors
In all simulations above, the bedrock component of the ice model is used, with non-local lithospheric deflection and local asthenospheric relaxation (PD09, PD12).That raises the possibility that bed elevations may have departed from modern observed values, and the inversion procedure could have erroneously compensated for biases in the bedrock model.Fortunately, in all simulations with the inverse method, bed elevations remain very close to those observed, with differences generally less than ∼ 20 m.As illustrated in Fig. D1, larger errors occur only in two isolated regions, northwesternmost Marie Byrd Land and northernmost Victoria Land, where the bed is ∼ 20 to 30 m too low; this is likely caused by isostatic depression under too thick ice over local mountains (same run shown in Fig. C1d-f).It seems likely that such small bedrock errors have had a very minor effect on the inverse results presented above.We tested this explicitly in one case, repeating the inverse run in Fig. C1d-f except with bedrock physics switched off so that bed elevations remain exactly at modern values.The results for surface elevations, basal coefficients and temperatures were essentially unchanged (not shown).Note, however, that bed elevation errors of this order might be significant for other purposes, for instance in comparing with relative sea level records.Also, note that the closeness of the agreement with observed modern bed elevations is somewhat fortuitous, because the model has not taken transient residuals from the last deglaciation into account (see Appendix E).

Appendix E Equilibrated versus transient modern state
In the inverse procedure, the model is forced with invariant modern climate, but its ice surface results are compared with modern observations.Any unequilibrated glacial isostatic adjustments (GIA) remaining from the last deglaciation are not accounted for in the model, and are implicitly assumed to be small.This mainly concerns ice mass inertia, lagged isostatic bed response, and ice temperatures through their effect on rheology and sliding.In future work we plan to combine the inverse method with transient runs through the last deglaciation, comparing in depth with relevant data (Sect.8; cf.Briggs et al., 2011;Whitehouse et al., 2012).For now, we can estimate the magnitude of this bias in the current results by comparing the model's equilibrated modern state with that at the end of a transient simulation through the last several 10 000s of yr.Results are shown in Fig. E1, where the transient simulation was run from 80 ka to the present, with paleoclimatic air temperatures, precipitation, oceanic melt rates and sea level parameterized as in PD12 (similarly to PD09).
As shown in Fig. E1a, equilibrated surface elevations on the East Antarctic central plateau are higher by ∼ 20 to 50 m than those in the transient run at 0 ka, presumably due to increasing snowfall rates over the last deglaciation.Conversely, surface elevations fall in many regions closer to the East Antarctic coast, generally coinciding with high bedrock topography and frozen basal temperatures.This may be due to continued warming of internal ice temperatures, lower ice viscosities and greater SIA shear flow.The marked increase in elevations over the Siple Coast is due to the model grounding line advancing slightly in the equilibrated run from its location at 0 ka in the transient run (seen in Fig. E1c vs. d), due to continuing bedrock rebound after the last deglaciation which causes grounding-line depths to shoal slightly.Changes in bed elevations (Fig. E1b) generally mirror the ice surface changes (smaller, with opposite sign, as expected due to isostatic relaxation), at least in East Antarctica.There is very little change in the frozen vs. melting pattern of basal temperatures (Fig. E1c-d).
We plan to assess all these effects as part of further transient modeling of the last deglaciation (Sect.8).Here, we note that the differences in ice surface elevations (Fig. E1a) are on the order of 20 to 50 m in most areas.These are on the same order as the residual large-scale elevation errors vs. observed achieved by the inverse method, and are much smaller than the many 100s of m elevation biases with simple C(x, y) distributions (Fig. 1a).Therefore, when the inversion method is extended in future work to properly account for transient effects, we anticipate that the deduced C(x, y) patterns will change only slightly in order to adjust for the elevation differences in Fig. E1a, and the overall C(x, y) results presented above will remain very much the same. in the inputs.The following 3 tests suggest that the deduced C(x, y) pattern is actually quite robust, with most regionalscale features insensitive to reasonable levels of uncertainty in the input fields.
The distribution of geothermal heat flux (GHF) under continental Antarctica is particularly uncertain.Commonly used datasets differ substantially from each other, and influence modeled basal temperatures on regional scales (Pattyn, 2010).Figure F1 shows the results of the inverse procedure using 3 different GHF distributions.Although these distributions cause significant regional differences in ice thicknesses and basal temperatures in non-inverse runs with this model (PD12), the inversion procedure accommodates the differences by quite small adjustments in the deduced C(x, y) distribution.With the low GHF values of Shapiro and Ritzwoller (2004) in the East Antarctic interior, C(x, y) values increase slightly in scattered areas, allowing more sliding to compensate for more basal freezing (Fig. F1g).However, on most regional and large scales, the patterns of deduced basal sliding coefficients are largely unaffected by the choice of GHF distribution (Fig. F1, bottom two rows).
Another uncertain input field is modern annual snow accumulation.The two datasets used here differ significantly around the ice-sheet margins (Fig. F2a and d).However, again, the resulting distributions of C(x, y) from the inverse procedure are very much the same.The only noticeable difference is in small patches in the East Antarctic Wilkes Land interior (Fig. F2g), compensating for slightly higher snow accumulation there in Arthern et al. (2006) compared to van de Berg et al. (2006).
Modern bedrock topography (LeBrocq et al., 2010) is used to calculate ice-free rebounded elevations, an input to the bed module of the ice-sheet model (PD12).In some regions, this topography is uncertain due to the sparsity of data lines (Lythe et al., 2001).Erroneous bed elevations could directly affect balance velocities, driving stresses, and potentially our deduced C(x, y) values.We crudely test this here by adding random Gaussian noise to the modern bedrock dataset, similarly to the random noise tests in Joughin et al. (2004).The noise is smoothed at small scales, retaining scales of 200 km or larger.With a root-mean-square noise amplitude of ∼ 200 m (Fig. F3c, second column), the effects on the deduced C(x, y) are limited to a few local areas such as the Siple Coast, and the large-scale patterns are scarcely affected (Fig. F3e and f).Increasing the noise amplitude to ∼ 400 m (Fig. F3g-j, third column) does begin to cause significant large-scale changes in C(x, y), for instance over the Gamburtsev Mountains.We conclude that the larger-scale C(x, y) patterns are robust to bed topographic errors up to a few hundred meters, and would be changed only by widespread errors of ∼ 400 m or more.

Appendix G Recovery Glacier Basin topography
Recently Le Brocq et al. (2008) suggested, based on observed ice surface curvatures, that much of the bed in the catchment of Recovery Glacier (a major system flowing into the Filchner Ice Shelf) may be much deeper than previously thought, up to ∼ 1500 m below sea level and layered with deformable sediments (Fig. G1).LeBrocq et al. (2011) found that the deeper topography improved modern surface elevations and velocities in their regional ice-sheet simulations (their Fig. 3).By applying the inversion procedure with both the standard and deeper topographies (both available in ALBMAP, Le Brocq et al., 2010), we can test if one or the other is more viable from the point of view of inversefitting to ice surface elevations.Figure G2 shows results us-ing the inverse method with sub-grid topographic influence, in a nested domain with 10-km resolution.
As might be expected, the deeper topography requires slower sliding velocities to compensate for the thicker ice (Fig. G2d and h), which is accomplished by the inversion procedure deducing somewhat less deformable sediment in the deeper Recovery and Slessor Glacier channels (Fig. G2b  and f).The surface elevation errors are generally slightly less with the "standard" topography, although both have errors of several hundred meters (of opposite sign) on the ridge between Slessor and Recovery Glaciers (Fig. G2a and e).Some of the elevation errors coincide with patches of frozen vs. melting bed (Fig. G2c and g), suggesting that the model's thermodynamics and parameterized effects of basal temperature and sliding might be at fault; however, much the same results are obtained with the inverse method that ignores   Brocq et al., 2008).Right column (i): observed surface ice velocity (Rignot et al., 2011), averaged here to 10-km model cells, m a −1 .

Fig. 1 .
Fig. 1.Non-inverse (left) vs. basic inverse (right) results.Top row: model minus observed surface ice elevation, meters.Middle row: basal sliding coefficients C(x, y), log 10 (m a −1 Pa −2 ).Bottom row: homologous basal temperatures T b (relative to pressure melt point), • C. Left column (a-c): using prescribed two-value C(x, y) according to whether rebounded ice-free bed elevations are above or below sea level.Right column (d-f): using basic inverse method.

Fig. 2 .
Fig. 2. Histograms of the magnitude of surface ice elevation errors| h s | (m) vs. number of grid cells over all grounded ice.Elevation error bins are in equal increments of log 10 (| h s |).Red: using twovalued basal sliding coefficient distribution as in Fig.1a.Green: using inverse method as in Fig.1d.Blue: using inverse method with sub-grid topographic influence as in Fig.3d.

Fig. 3 .
Fig. 3. Improved results with sub-grid basal topographic influence on sliding.Top row: model minus observed surface ice elevation, meters.Middle row: basal sliding coefficients C(x, y), log 10 (m a −1 Pa −2 ).Bottom row: homologous basal temperatures T b , • C. Left column (a-c):using basic inverse method, same as in Fig.1d-f.Middle column (d-f): using inverse method with sub-grid topographic influence (s.a.).Right column (g-i): non-inverse run using prescribed C(x, y) from (e) with sub-grid topographic influence, and freely varying grounding lines and ice sheets.

2
. model errors in basal mass balance (melting, refreezing) are negligible, 3. the model ice thicknesses are correct, 4. the model's split between surface and depth-average velocities is correct, and 5. both model and real modern Antarctic ice sheets are close to equilibrium (Appendix E).

Fig. 4 .
Fig. 4. Time series of mean absolute surface elevation error (model vs. observed | h s |, m) in inverse runs with sub-grid topographic influence, for various combinations of t inv (time interval between adjustments, yr) and h inv (scaling constant, m) in Eq. (5).

Fig. 5 .
Fig. 5. (a) Observed surface ice velocities(Rignot et al., 2011), averaged here to 20-km model cells, m a −1 .(b) Model surface ice velocities, m a −1 , using inverse method with sub-grid topographic effect, 20-km resolution (as in Fig.C1d-f).(c) Model minus observed velocities, log 10 (m a −1 ), i.e., log 10 (v model /v observed ).Very slow velocities are ignored, i.e., if v model or v observed is less than 2 m a −1 , it is reset to 2 m a −1 for this plot.(d) Scatter plot of observed vs. model velocities, log 10 (m a −1 ), for each 20-km grid cell with grounded ice.The same figure appears in PD12.

Fig. A1 .
Fig. A1.Poor results using inverse method with no basaltemperature effect on sliding.Top row: model minus observed surface ice elevation, meters.Middle row: basal sliding coefficients C(x, y), log 10 (m a −1 Pa −2 ).Bottom row: homologous basal temperatures T b , • C. Left column (a-c): using basic inverse method (no sub-grid topographic influence) and with no effect of basal temperature T b on sliding.Right column (d-f): non-inverse run using prescribed C(x, y), nearest-neighbor filled from (b).

Fig. B1 .
Fig. B1.Crude bounding of enhancement factor for internal flow (E).Top row: model minus observed surface ice elevation, meters.Middle row: basal sliding coefficients C(x, y), log 10 (m a −1 Pa −2 ).Bottom row: homologous basal temperatures T b , • C. Each column shows results using basic inverse method (no sub-grid topographic influence) for a different value of the internal-flow enhancement factor E, ranging from E = 0.1 to E = 8.

Fig. E1 .
Fig. E1.Equilibrated vs. transient model simulations at 0 ka.(a) Difference in ice surface elevations (meters) between an equilibrated run with invariant modern forcing minus a snapshot at 0 ka from a transient run through the last 80 kyrs.Both runs are noninverse, using prescribed inverse-derived basal coefficients C(x, y) (Fig. 3e) and freely varying grounding lines and ice shelves as in Fig. 3g-i.(b) As (a), except difference in bed elevations.(c) Homologous basal temperature ( • C) for transient run at 0 ka.(d) As (c), except for equilibrated run.

Fig. G2 .
Fig. G2.Results with standard (left) vs. deeper (right) bed topography in the Recovery Glacier basin.Top row: model minus observed surface ice elevation, meters.Second row: basal sliding coefficients C(x, y), log 10 (m a −1 Pa −2 ).Third row: homologous basal temperatures T b , • C. Bottom row: surface ice velocity, m a −1 .Left column(a-d): using inverse method with sub-grid topographic influence, limited-domain nested run at 10-km resolution with lateral boundary conditions from corresponding continental run.Middle column (e-h): as (a-d), except using alternate bed topography with deeper elevations in Recovery Glacier drainage basin (LeBrocq et al., 2008).Right column (i): observed surface ice velocity(Rignot et al., 2011), averaged here to 10-km model cells, m a −1 .