Geoscientific Model Development Description of a hybrid ice sheet-shelf model , and application to Antarctica

The formulation of a 3-D ice sheet-shelf model is described. The model is designed for long-term continentalscale applications, and has been used mostly in paleoclimatic studies. It uses a hybrid combination of the scaled shallow ice and shallow shelf approximations for ice flow. Floating ice shelves and grounding-line migration are included, with parameterized ice fluxes at grounding lines that allows relatively coarse resolutions to be used. All significant components and parameterizations of the model are described in some detail. Basic results for modern Antarctica are compared with observations, and simulations over the last 5 million years are compared with previously published results. The sensitivity of ice volumes during the last deglaciation to basal sliding coefficients is discussed.


Introduction
This paper describes the formulation of a combined ice sheetshelf model, some aspects of which have been included in earlier papers (Pollard andDeConto, 2007, 2009; henceforth PD07, PD09), but many have not.Here, a full model description is presented, including recently added features that are being used in current work (Pollard and DeConto, 2012;henceforth PD12).
Early numerical 3-D ice sheet models used the shallow ice approximation (SIA), i.e., scaled dynamical equations appropriate for large-scale grounded ice flow dominated by vertical shear ("∂u/∂z") and basal stress locally balancing the gravitational (surface-slope) driving stress (e.g., Andrews and Mahaffy, 1976;Budd and Smith, 1979).Later, the need to include floating ice shelves, fast flowing ice streams (with very low basal drag) and grounding-line migration emerged, for instance to model marine collapses of the West Antarctic Ice Sheet.For shelves and streams, flow is dominated by horizontal stretching ("∂u/∂x") and a different set of scaled equations is appropriate using the shallow shelf approximation (SSA, also called shelfy stream) (Morland, 1982;MacAyeal, 1989MacAyeal, , 1996)).This was first attempted in 3-D models in the late 1990's, by applying either the SIA or SSA equations in different specified regions with matching at the boundaries (Hulbe and MacAyeal, 1999;Huybrechts and de Wolde, 1999;Huybrechts, 2002;Ritz et al., 2001;cf. Budd et al., 1994).
However, ice-stream flow can consist of both vertical shear and horizontal stretching, and its boundaries are not always amenable to simple parameterization; also, more rigorous treatment of the grounding zone is needed for accurate simulations of grounding-line migration (Schoof, 2007).Models with more rigorous (i.e., less scaled) higher-order or full-Stokes formulations of the dynamical equations are available (e.g., Pattyn, 2002;Seddik et al., 2012), which capture both modes of flow, but are computationally expensive and are not currently feasible for long-term continental-scale applications.The problem has been addressed by a number of hybrid models, which use heuristic superpositions of the depth-integrated SIA and SSA equations (Bueler and Brown, 2009;Winkelmann et al., 2011;Goldberg, 2011;formalism in Schoof and Hindmarsh, 2010; different approaches in Marshall and Clarke, 1997;Hubbard, 2006;early steps in Alley and Whillans, 1984;van der Veen, 1985).Although these models are not rigorous, their performance can be tested against higher-order/full-Stokes models in simple scenarios as noted below, and are feasible for long-term large-scale applications.
The model described here is a hybrid model, most similar to Winkelmann et al. (2011) and Goldberg (2011).Our dynamical equations fall into type L1L2 of Hindmarsh's (2004) categories.An additional measure is taken to improve simulation of grounding zones, where the Schoof (2007) parameterization is imposed as a condition on ice flux across the grounding line.This enables grounding-line migration to be simulated reasonably accurately without the need for much higher resolution (Schoof, 2007;Gladstone et al., 2010aGladstone et al., , 2012a;;Pattyn et al., 2012a).
The model also includes standard equations for the evolution of ice thickness, internal ice temperatures, and the bedrock response under the ice load.An optional coupling with a sediment model, with explicit quarrying/abrasion, transport and deposition of deformable sediment under the ice, is fully described in Pollard andDeConto (2003, 2007) and is not covered here.There is no explicit basal hydrologic component in the current model.
For reference, new features added to the model since PD09 and described below are listed here: -new parameterization of oceanic melt at base of floating ice; -calving parameterization at floating ice edge; -sub-grid fractional area of ice vs. ocean in cells at floating ice edge; -oceanic melting at vertical ice faces; -parameterization of shelf drag by sub-grid bathymetric pinning points; -modified sub-grid application of Schoof (2007) grounding-line condition; -optional simplifications in the combined SIA-SSA dynamics; -adaptive reduction of model time step to avoid numerical instability; -distribution of basal sliding coefficients deduced by a simple inverse method, described in PD12, and with the resulting pattern used here.
Two other features, not used in the applications below, will be described in future papers: -sub-grid ice surface elevation interpolation and fractional area for calculation of surface mass balance at terrestrial ice margins (cf.van den Berg et al., 2006); -improved numerics for nesting model capability in higher-resolution limited domains, with lateral boundary conditions from a previous continental run.
The bulk of this paper (Sects.2.1 to 2.13) contains the model description, followed by an account of input datasets and climate forcing in Sect.3. Section 4 presents results for modern Antarctica, where simulations at different resolutions are compared with observations.Section 5 presents paleoclimatic simulations of the last 5 Myrs, repeating those in PD09 with the new model version, and briefly discusses issues concerning the last deglaciation.

Model description
The model consists of diagnostic equations for ice velocities, and 3 prognostic equations for the temporal evolution of ice thickness, ice temperature, and bedrock deformation below the ice.Prescribed boundary fields are equilibrium bedrock topography and corresponding loading (modern rebounded ice-free state), unfrozen basal sliding coefficients, geothermal heat flux, and sea level.Monthly mean surface air temperatures and precipitation are either parameterized or provided from a climate model, in order to calculate annual surface mass balance and ice surface temperature (there is no seasonal cycle in the ice model itself).Sub-ice oceanic melting and shelf-edge calving are parameterized for floating ice shelves.A list of model symbols is provided in Table 1.

Horizontal and vertical grids
The ice sheet-shelf model uses a finite-difference Arakawa-C grid (e.g., Rommelaere and Ritz, 1996), where horizontal velocities (u, v) are calculated on separate grids staggered by half a grid box relative to ice thickness (h), as shown in Fig. 1.The model code contains metric terms appropriate for Cartesian, Polar Stereographic, and Spherical Polar (longitude-latitude) grids, and also for flow lines with one horizontal dimension.Note however that for longitudelatitude grids, a rigorous derivation of the SSA equations introduces some spherical metric terms not in the current code, Total ice velocities, x and y directions, across grounding line (m a −1 ) q g Total ice flux, x or y direction, across grounding line (m 2 a −1 ) Basal ice homologous temperature ( Ice flow enhancement factor (1 for SIA, 0.3 for SSA) which would need to be modified in order to properly treat global-scale floating ice (Tziperman et al., 2012).
The ice model uses a vertical coordinate z running from 0 at the ice surface to 1 at the ice base: where h s is ice surface elevation and h is ice thickness.The vertical grid has 10 uneven layers, more closely spaced near the top and bottom.Ice temperatures and horizontal velocities are defined at the mid point of each layer.

Ice velocities
A combination of the scaled equations for vertical shearing ("∂u/∂z", shallow ice approximation, SIA) and horizon-tal longitudinal stretching ("∂u/∂x", shallow shelf approximation, SSA) is used, similarly to several previous hybrid models mentioned in Sect. 1.The combination is heuristic because neither scaling is accurate where both shearing and stretching are significant (streaming and grounding zones).Nevertheless, with the additional imposition of Schoof's (2007) grounding-line flux condition described below, reasonable results are obtained in idealized intercomparison tests (Pattyn et al., 2012a, b).The numerics are simple enough to allow long-term O(10 7 ) year continental-scale runs.Recent modeling progress using full-Stokes or higherorder flow equations on fine or adaptive grids rigorously include these processes (e.g., Morlighem et al., 2010;Seddik et al., 2012), but require considerably more computer time, and for now are limited to shorter time or smaller spatial scales.

P
Annual mean precipitation rate (m a −1 ice equivalent) C Basal sliding coefficient between bed and ice (m a −1 Pa −2 ) C(x, y) Basal sliding coefficient for unfrozen beds (m a −1 Pa −2 ) C froz Basal sliding coefficient for no flow (10 −20 m a −1 Pa −2 ) m Basal sliding exponent (2) T r Threshold temperature in basal sliding (−3 Latent heat of fusion (0.335 × 10 6 J kg −1 ) q Bed load (Pa) As described in PD07 and PD09, the SIA and SSA equations are combined as follows: 1.In the expressions for effective viscosities, SIA's ∂u/∂z shear-softening terms are included in the viscosity for SSA, and SSA's ∂u/∂x terms are included in the viscosity for SIA.
2. The SSA equations solve for depth-averaged total velocity.But in the SSA basal stress terms, a distinction is made between depth-averaged and basal velocity, with the difference being the vertical mean of the SIA shear flow.
3. The driving stress in the SIA equations is reduced by the gradient of the longitudinal stress from the SSA equations acting on the column above each level.
These steps require an iteration between the SSA and SIA solutions, as described below.Goldberg (2011) takes essentially the same steps, and discusses the relationship with Schoof and Hindmarsh (2010).Cartesian coordinates are used in the equations here, although metric terms are included in the model code to handle other grids (but see Sect.2.1).All symbols are listed in Table 1.The following presentation is very similar to PD07 Appendix A (noting an erroneous factor of 2 in their Eq.A2a, b).
Writing Cartesian horizontal ice velocities as u(x,y,z) and v(x, y, z), define the basal ice velocity u b (x, y) = u(x, y, z b ), and the internal shearing ice velocity u i (x, y, z) = u − u b , so that u i (x, y, z b ) = 0. Denoting vertical averages through the ice column with a bar, then ū = u b + ūi and similarly for v, v b and vi .The SIA-like internal shear equations for u i (x, y, z) and v i (x, y, z) are where σ ij are deviatoric stresses given below (Cuffey and Paterson, 2010).The SSA-like horizontal stretching equations for ū(x, y) and v(x, y) are Equation (2a, b) and their horizontal boundary conditions for unconfined ice shelves are derived for instance in Morland (1982) and MacAyeal (1996).In Eqs.(2a, b), u b = ū − ūi and v b = v − vi where ūi and vi are obtained from vertical integrations of (1a, b) (e.g., Ritz et al., 1997).
In the zero-order shallow ice approximation, the vertical shear stress (σ xz , σ yz ) in Eqs.(1a, b) would be balanced only by the hydrostatic driving force −ρ i g(h s − z) (∂h s /∂x, ∂h s /∂y) acting on the ice column above level z.
Here, horizontal stretching forces are included in this force balance (Hubbard, 1999(Hubbard, , 2006;;Marshall et al., 2005), so that where LHS x and LHS y are the left-hand sides of Eqs.(2a) and (2b), respectively.Because horizontal stretching forces are taken to be vertically uniform and the terms in Eq. ( 2) are forces on the entire ice thickness, their effect on the ice column above level z is scaled by (h s − z)/h in Eq. (3).Inclusion of the strain softening terms in Eqs. ( 1) and (2) due to each other's flow (i.e., σ xx , σ yy , σ xy in Eq. ( 1), ∂u i /∂z and ∂v i /∂z in µ in Eq. ( 2)), requires manipulation of the constitutive relation for ice rheology.In Eq. (2), and A = ∫ A dz/h is the vertical mean of the Arrhenius temperature-dependent coefficient in the constitutive relation where εij are strain rates, σ ij are deviatoric stresses, and ε and σ are the second invariants of their respective tensors.The latter are defined by ε2 ≡ The corresponding expression for σ 2 is used in Eq. ( 1), and the purely horizontal components are obtained in our numerical procedure from The basal sliding relation used on the right-hand sides of Eqs.(2a) and (2b) for grounded ice is ũb , where τb is basal stress.Where ice is grounded, i.e., where ρ w (S −h b ) < ρ i h or the ocean has no access (held back by intervening thicker ice or higher land), then f g =1 in the sliding terms, and the ice surface elevation h s = h + h b .Where ice is floating, i.e., ρ w (S − h b ) > ρ i h and the ocean has access, then At each time step, an outer iteration is performed that solves for SSA and SIA velocities, updates ice thicknesses for half of the time step, and re-solves the velocities using the new ice thicknesses, i.e., a second-order Runge-Kutta method.In the solution of Eq. (2) for SSA velocities, a standard (Picard) inner iteration is performed to account for the non-linear dependence of µ and basal sliding on the velocities in Eqs. ( 2), ( 4) and ( 6).The outer iteration converges naturally to the appropriate scaling of SSA vs. SIA flow, depending on the magnitude of the basal sliding coefficient.Usually the flow is either almost all vertical shear, with basal drag balancing the driving stress and with negligible stretching, or is almost all longitudinal stretching which balances the driving stress, with small or no basal drag and negligible internal shear.For a fairly narrow range of sliding coefficients, significant amounts of both flow types co-exist.
In each pass of the outer iteration, the SSA Eqs.(2a, b) are solved first, using a sparse matrix method, or optionally, successive over-relaxation (SOR) (or a tridiagonal matrix solution for 1-D flow-line problems).Then the ice-thickness advection equation (Sect.2.6) is time-stepped accounting for both SSA and SIA flow.Ice advection due to SIA is performed time implicitly, with the vertically averaged SIA flow given from Eqs. ( 1) and ( 3) and using time-implicit linearized Newton-Raphson contributions from all h and ∂h s /∂x terms (as in earlier SIA-only model versions; DeConto and Pollard, 2003).Centered ice thicknesses are used for the SIA advection, whereas the time-explicit SSA advection uses upstream ice thicknesses for stability (in Eq. 14 below).An Alternating Direction Implicit (ADI) scheme is used for x vs. y directions (Mahaffy, 1976).A CFL-based maximum speed limit on horizontal velocities can be imposed for stability.No ice advection is allowed out of grid cells with sub-grid areal fraction f e < 1 (which occurs only for cells at the edge of floating ice shelves, see Sect.2.9).
CPU time in the model is dominated by the Sparse-Matrix (or SOR) solutions of the SSA equations.As described in PD09, a considerable reduction in CPU time can be achieved by restricting the full SSA+SIA iterative procedure to grid points with mid-to-high values of the basal sliding coefficient, C(x,y) ≥ 10 −8 m a −1 Pa −2 (see PD12).This range includes all fast streaming regions underlain with deformable sediment (∼10 −5 ).For lower C(x, y) values < 10 −8 (including hard bedrock, ∼10 −10 ), the full procedure yields virtually 100 % shearing (SIA) flow anyway.At the latter points, advection by internal shear deformation ( ūi , vi ), and basal sliding (u b , v b ) are both modeled by standard SIA dynamics.At full SSA+SIA points with C(x,y) ≥ 10 −8 , advection by internal shear deformation, basal and horizontal stretching are all included in the coupled solution of Eqs. ( 1) and ( 2).Tests show that results are essentially unchanged from those with the full SSA+SIA iteration performed at all locations.
In intermediate model versions, some simplifications were tried in the coupling dynamics such as neglecting the strain softening cross-terms in Eqs. ( 1) and ( 6), which reduced CPU time modestly with only slight effects on the results.Some of these simplifications were used for the figures shown below; however, the most complete and current model version is described above.

Grounding-line flux condition
Flow-line tests with hybrid or higher-order models show that in order to capture grounding-line migration accurately, it is necessary either to resolve the grounding-zone boundary layer at very fine resolution (Schoof, 2007;Gladstone et al., 2010aGladstone et al., , 2012a;;Pattyn et al., 2012a), or to apply an analytic constraint on the flux across the grounding line.The latter approach is used here, with flux q g across model grounding lines parameterized as in Schoof (2007, his Eq. 29): This yields the vertically averaged velocity u g = q g /h g where h g is the ice thickness at the grounding line.The middle term in Eq. ( 8) accounts for back stress at the grounding line due to buttressing by downstream islands, pinning points or side-shear, where τ xx is the longitudinal stress just downstream of the grounding line, calculated from the viscosity and strains in a preliminary SSA solution with no Schoof constraints.The free stress τ f is the same quantity in the absence of any buttressing, given by 0.5 ρ i gh g (1−ρ i /ρ w ) (cf.Goldberg et al., 2009;Gagliardini et al., 2010).Ā is the depth-averaged ice rheological coefficient and n is the Glen-Law exponent, C s is Schoof's (2007) basal sliding coefficient and m s the basal sliding exponent, corresponding here to C −1/m and 1/m, due to the reversed form of the basal sliding law.ρ i and ρ w are densities of ice and ocean water respectively, and g is the gravitational acceleration.h g is interpolated in space by first estimating the sub-grid position of the grounding line between the two surrounding floating and grounded h-grid points.This is done by linearly interpolating height above flotation between those two points to where it is zero, linearly interpolating bedrock elevation to that location, and then simply computing the flotation thickness of ice for that bedrock elevation and current sea level (equivalent to LI in Gladstone et al., 2010b).The velocity u g is calculated at the grounding-line points on the u-grid, i.e., those with floating ice in one adjacent (left or right) h-grid box and grounded ice in the other (and similarly for v g on the v-grid).These velocities are imposed as an internal boundary condition for the flow equations, in effect overriding the large-scale velocity solution at the grounding line.This procedure only considers one-dimensional dynamics perpendicular to the grounding line, as in the 1-D flowline analysis in Schoof (2007).It works naturally with the staggered C-grid (Sect.2.7), where the grounding line is a continuous series of perpendicular segments of u-direction or v-direction interfaces between h-grid boxes, and u g (v g ) velocities flow across interfaces running through u-grid (vgrid) points.Equations (8-9) and the discussion in this section applies equally to the y direction, with v g and τ yy instead of u g and τ xx .Note that spatial gradients of quantities parallel to the grounding line, which are not included in Schoof's (2007) flow-line derivation of Eq. ( 8), are neglected here (Katz and Worster, 2010;Gudmundsson et al., 2012;Pattyn et al., 2012b).
We have tested this method of solution in many idealized 1-D flow-line tests, similar to those in Schoof (2007).Our goal was to achieve the same grounding-migration results using a coarse grid (∼10 to 40 km) with those using very finegrids (∼ 0.1 km).For coarse grids, we find that it is necessary to impose Eq. ( 8) as a grounding-line boundary condition.Also for coarse grids we find that an additional rule is necessary, because the outer-solution structure of the grounding zone is not fully captured by the grid: 11 flotation thickness of ice for that bedrock elevation and current sea level (equivalent to LI in Gladstone et al., 2010b).
The velocity ug is calculated at the grounding-line points on the u-grid, i.e., those with floating ice in one adjacent (left or right) h-grid box and grounded ice in the other (and similarly for vg on the v-grid).These velocities are imposed as an internal boundary condition for the flow equations, in effect overriding the large-scale velocity solution at the grounding line.This procedure only considers one-dimensional dynamics perpendicular to the grounding line, as in the 1-D flowline analysis in Schoof (2007).It works naturally with the staggered Cgrid (section 2.7), where the grounding "line" is a continuous series of perpendicular segments of u-direction or v-direction interfaces between h-grid boxes, and ug (vg) velocities flow across interfaces running through u-grid (v-grid) points.Eqs.(8-9) and the discussion in this section applies equally to the y direction, with vg and τyy instead of ug and τxx.Note that spatial gradients of quantities parallel to the grounding line, which are not included in Schoof (2007)'s flowline derivation of Eq. ( 8), are neglected here (Katz and Worster, 2010;Gudmundsson et al., 2012;Pattyn et al., 2012b).
We have tested this method of solution in many idealized 1-D flowline tests, similar to those in Schoof (2007).Our goal was to achieve the same grounding-migration results using a coarse grid (~10 to 40 km) with those using very fine-grids (~0.1 km).For coarse grids, we find that it is necessary to impose (8) as a grounding-line boundary condition.Also for coarse grids we find that an additional rule is necessary, because the outer-solution structure of the grounding zone is not fully captured by the grid: If the flux qg from Eq. ( 8) is greater than the large-scale shelf-equation's flux qm at the grounding line, then ug (= qg/hg) is imposed exactly at the u-grid grounding-line point; conversely if qg < qm, then ug is imposed one u-grid box downstream of the grounding-line point.The former is usually associated with grounding-line retreat, and the latter usually with grounding-line advance.
(9) (9) When converting the grounding-line flux q g from Eq. ( 8) to a velocity (u g ), it is important to divide by the ice thickness (called h g above) that will effectively be used at the relevant point in Eq. ( 9) in the finite-difference numerics of the ice advection equation.Then the model's flux at that point will be exactly that from Eq. ( 8).In simple equilibrated flow-line tests, this means that the model flux equals the net surface mass balance upstream from the grounding line, an important property of analytic solutions.This yields good agreement with analytic solutions including hysteresis in MISMIP-like tests, using grid sizes of ∼ 5 to several 10s km (Pollard and DeConto, 2011;Docquier et al., 2011;Pattyn et al., 2012a).The agreement can be made almost exact by adjusting the flux q g for the increment in surface mass balance between the actual grounding line and the point where Eq. ( 9) is applied, as illustrated in Fig. 2. The analytic solutions in turn agree well with full-Stokes model results, at least in steadystate non-transient situations (Drouet et al., 2011;Pattyn et al., 2012a).
In efforts to minimize single-cell dithering in some idealized tests, i.e., flipping back and forth between upstream and downstream points in Eq. ( 9), two further measures were taken: i.An initial SSA solution is done at each time step, without any imposed flux from Eq. ( 8), to calculate the large-scale flux that is compared to the imposed flux in Eq. ( 9).Previously the large-scale flux was estimated by local finite differences.
ii.Values of the imposed velocities from Eq. ( 8) are calculated for both upstream and downstream points of the grounding line, and these are imposed in the flow equations with weights between 0 and 1 depending on how much (and with what sign) the large-scale flux differs from the imposed flux as in Eq. ( 8).
These measures had little effect on the dithering in flow-line tests, but fortunately no associated degradation of large-scale results has been detected.

Basal sliding coefficients
Basal sliding is treated similarly to PD09 by a standard drag law (Cuffey and Paterson, 2010;Pattyn, 2010;Le Brocq et al., 2011) where ũb is basal sliding velocity, τ b is basal stress, and m = 2 as in Sect.2.2.As described in PD12, the sliding coefficient C' depends on homologous basal temperature, implicitly representing basal hydrology: ) where r = max 0, min 1, (T b + 3) 3 where C(x, y) is the full sliding coefficient, and C froz = 10 −20 ma −1 Pa −2 (which is small enough to prevent any discernible sliding, but is not exactly zero to avoid divide-byzero exceptions in the numerics).T b ( • C) is the homologous basal temperature, i.e., relative to the pressure melt point T m = −.000866h where h is ice thickness (m).There is no sliding below the threshold homologous temperature (−3 • C), ramping up linearly to full sliding at the melt point.
C(x, y) is a specified basal sliding coefficient representing intrinsic bed properties.In PD09 it was twovalued, depending on whether the modern rebounded Antarctic bedrock is above or below sea level: if above, C(x, y)=10 −10 m a −1 Pa −2 representing hard bedrock (mainly EAIS), and if below, C(x, y) = 10 −6 m a −1 Pa −2 representing deformable sediment (mainly WAIS) (e.g., Studinger et al., 2001) shown in Fig. 3a.In PD12, a simple inverse method is used that attempts to deduce the real distribution of C(x, y) under modern Antarctica, constrained to the range 10 −10 to 10 −5 m a −1 Pa −2 .
In PD12, modern Antarctic results are further improved by adding a dependence on sub-grid bedrock relief, that allows more sliding across major mountain ranges, presumably in deep and warmer valley troughs not resolved by the model grid.Without this addition, basal ice is often completely frozen over mountain ranges, and insufficient crossrange flow causes surface elevations to be too high (PD12).We attempt to parameterize this sub-grid process by modifying the width of the basal-temperature ramp in Eq. ( 11), replacing it 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.This quantity was also used by Marshall et al. (1996) in another context.h eq b is the ice-free, isostatically rebounded, 9-pointsmoothed bed elevation on the model grid, used to mimic SA in data-sparse regions (PD12).The values of the constants are discussed in PD12.Whitehouse et al. (2012) apply a similar increase in sliding coefficient over mountainous terrain, for much the same reasons.Equation ( 12) and the associated inverse-derived C(x, y) distribution (Fig. 3b) are used in the simulations below.
For grid points where the full SSA+SIA iteration is performed (Sect.2.2), C and u b enter in the right-hand side of the SSA in Eq. ( 2), where Eq. ( 10) is inverted to give τ b as a function of u b , and u b is treated time explicitly in the stepping of the ice thickness as in Eq. ( 14).For points where just the SIA equation is used as discussed in Sects.2.2, Eq. ( 10) is treated time implicitly, with τ b equal to the driving stress (ρ i g h ∂h s /∂x), and with linearized Newton-Raphson contributions from Eq. ( 10) in the same way as for internal shearing in Eq. (1).

Sub-grid pinning points
Under the major ice shelves, there may be sub-grid pinning points due to small bathymetric rises scraping the bottom of the ice, especially near the grounding line, that are unresolved by the model grid.This is parameterized simply in terms of the standard deviation of observed bathymetry within each model cell.The fractional area f g of ice in contact with sub-grid bathymetric high spots is where h w is the thickness of the ocean column between the cell-mean bedrock and ice base, and s dev is the standard deviation of the observed bed elevations (ALBMAP, 5 km, Le Brocq et al., 2010) within the cell.For 20 to 40 km grids, s dev is typically smaller than ∼50 m under the Ross and much of the Weddell and Amery ice shelves, but up to to a few 100s m in isolated patches of the Weddell, Lambert, and much of Pine Island Bay.f g here is identical to the f g in Eqs.(2a, b), and modifies the basal stress for the cell.Instead of no drag (f g = 0, freely floating ice), the value from Eq. ( 13) is used, increasing the basal stress to f g times the amount for 100 % basal contact.
In effect, this augments the overall drag on the ice shelf in addition to side drag, which is transmitted upstream within the SSA equations, increasing buttressing and reducing ice flux across the grounding line, i.e., making τ xx less positive and reducing q g in Eq. ( 8).The extent and importance of small-scale pinning is somewhat speculative, and deserves more study, observationally by examination of surface fea-tures or improved bathymetry (Horgan and Anandakrishnan, 2006;Fricker et al., 2009;Hulbe et al., 2010;Timmermann et al., 2010), and by modeling such as Favier et al. (2012).The time stepping of the ice thickness equation is done as part of the iterative solution of ice velocities as described in Sect.2.2.The treatments of the various local ice gains or losses (SMB, etc) are described in later sections.

Ice temperature and rheology
The prognostic equation for internal ice temperatures T (x, y, z , t) is where z = (h s − z)/ h, k i = 2.1 × 365 × 86 400 J a −1 m −1 K −1 is ice thermal conductivity, and Q i is internal shear heating due to both SIA and SSA deformation.Only vertical heat diffusion is included; horizontal heat diffusion is assumed negligible on scaling grounds.Note that the vertical coordinate z is dimensionless, and Eq. ( 15) has been transformed to this coordinate system (Huybrechts and Oerlemans, 1988;Ritz et al., 1997).The transformed vertical velocity w = dz /dt; numerical calculation of w uses the technique in Ritz et al. (1997) (whose w t is our w h).
Horizontal velocities u, v are the sum of internal (∼SIA) shear and the basal velocity.The large-scale advective terms (−u∂T /∂x -v∂T /∂yw ∂T /∂z ) are calculated timeexplicitly, using upstream parabolic interpolation for T (Farrow and Stevens, 1995).
The upper boundary condition is T(x,y,0,t) = surface ice temperature, deduced from surface air temperatures (Sect.3).For grounded ice, the lower boundary condition at the ice base is that the vertical conductive flux (k i / h) ∂T /∂z at  15) is time-stepped with the vertical diffusive terms and boundary conditions treated time implicitly, which involves a standard tridiagonal solution versus z for each ice column.To avoid numerical instability, very small ice thicknesses (< 1 m) are treated as a thin film with zero heat capacity, but still with latent heat and melting if its temperature would otherwise exceed the pressure melt point.
Surface melting, refreezing and locally mobile liquid are calculated along with the surface mass balance (Sect.3).Any locally mobile liquid (rain, snow melt and ice melt, minus refreezing) is assumed to immediately percolate downwards into the local vertical ice column, exchanging its latent heat with the sensible heat of the next lowest layer, i.e., if the layer is below freezing, then some (all) of the percolating liquid freezes, raising the layer temperature to (towards) the pressure melt point (and adding to the layer thickness).If the melt point is reached, the remaining water percolates down to the next layer, and so on.If any liquid reaches the base, it is added to any ice melt at the base itself, and is simply recorded as mass lost from the model (there is no basal hydrologic component).
The model includes vertical heat diffusion and storage in bedrock below the ice, heated from below by a specified geothermal heat flux.Nominally, and in all simulations shown below, its effect is minimized by using a very thin (30 m) single layer, so that the geothermal heat flux is essentially applied to the base of the ice.In other applications, it is typically ∼ 2 km thick with 6 unequally spaced layers (cf.Ritz et al., 1997).Physical and thermal properties of bedrock are given in Table 1.
In the ice dynamics (Sects.2.2 and 2.3), the ice rheological coefficient A and its dependence on temperature are specified as in Huybrechts (1998): where T is the homologous ice temperature T − T m , where T m = −.000866z is the pressure melting point ( • C) and z is depth (m) below ice surface.Units of A are a −1 Pa −3 corresponding to n=3 in Eqs.(1) to ( 7).The enhancement factor E is set to 1 for SIA flow in Eq. (1) (see PD12), and to 0.3 for SSA flow (Eqs. 2 and 8).The ratio of enhancement factors represents differences in fabric anisotropy between grounded and shelf ice (Ma et al., 2010); it is similar to but somewhat smaller than their suggested range of 5:1 to 10:1.

Sub-ice-shelf oceanic melting
The simulation of oceanic melting at the base of Antarctic ice shelves is challenging, involving incursions of Circumpolar Deep Water (CDW) or High Salinity Shelf Water (HSSW) and other mechanisms that differ from basin to basin (e.g., Nicholls et al., 2009;Walker et al, 2009;Jenkins et al., 2010;Pritchard et al., 2012).Coupling with ice sheet models will ultimately require high-resolution 3-D regional ocean modeling (e.g., Dinniman et al., 2011;Hellmer et al., 2012), especially for paleo and future scenarios.For now, we use simple parameterizations that attempt to provide (i) the basic modern spatial distribution, and (ii) paleoclimatic variations that yield results in accord with geologic data.
In PD09, the parameterization of modern oceanic melt rates was somewhat ad hoc, based on subtended arcs to open ocean.Our current parameterization described below follows Martin et al. (2011) for the PISM-PIK model.A new parameterization based on Olbers and Hellmer's (2010) more physical cavity-box model (Gladstone et al., 2012b;Winkelmann et al., 2012) is under development.
Similarly to Martin et al. (2011), the oceanic melt rate at the base of floating ice (m a −1 ), OMB in Eq. ( 14), is given by where T o is the specified ocean water temperature, and T f = .0939-.057×34.5 -.000764 z ( • C) is the ocean freezing point at ice-base depth z (m) (Beckmann and Goose, 2003;cf. Jenkins and Bombosch, 1995).The transfer factor and K is an additional O(1) basin-dependent factor given below.Because the freezing point T f decreases with depth, the dependence on T o − T f means that melt rates tend to be higher at the grounding line as deduced from observations.Unlike Martin et al. (2011), the dependence on temperature difference T o −T f is quadratic (Holland et al., 2008).
Here, the ocean temperature T o is specified differently for various Antarctic sectors, based on observations but mainly aiming to produce realistic modern ice-shelf extents and grounding-line positions.The 4 sectors are delineated by crude latitude and longitude ranges, as follows (with latitudes in • N, longitudes in • E, temperatures in • C, and depths in meters), and also shown in Fig. 4a.
-Amundsen and Bellingshausen Seas, and Western Peninsula: T o − T f depends on depth z, based loosely on profiles in the outer Pine Island Bay with an upper layer of colder fresher water (Jacobs et al., 2011), which may be important for the survival of smaller shelves with shallow grounding lines: T o − T f = 0.5 for z < 170, 3.5 for z > 680, joined linearly from 170 to 680 m.T o − T f and K are as for the Amundsen/Bellinghausen/W.Peninsula sector, even though ocean profile data in Prydz Bay for instance do not indicate a distinct upper layer as clearly as for Pine Island Bay (Smith et al., 1984).
-Ross embayment: where This has the effect of reducing ocean melting for regions mostly surrounded by land.It is found to be necessary in long-term paleo runs (Sect.5 below) to allow WAIS to regrow after a collapse of all marine ice.After a collapse, the surviving small terrestrial ice caps on Western Antarctic islands must first form thin ice shelves that grow over the interior seaway, coalesce, thicken and become buttressed so as to allow grounding lines to advance out from the islands.Equation ( 18) can be justified by arguing that interior seaways mosly surrounded by land were more protected from warm water intrusions than the modern coast and embayments.This hypothesis should be tested by regional ocean modeling of the environment following a major WAIS collapse.Equation ( 18) is not applied to East Antarctica for the ad hoc reason that ocean melting from Eq. ( 17) needs to penetrate into the Lambert Graben in order to produce reasonable modern grounding line and shelf extents there.A similar parameterization to Eq. ( 18) is also used to restrict calving (Sect.2.10).
The above yields the distribution of modern ocean melt rates, shown in Fig. 4b.For paleoclimatic applications, longterm climate variations are parameterized much as in PD09, based on a single weighting parameter w c set proportional to deep-sea core δ 18 O, plus a slight influence of austral summer insolation: where S is eustatic sea level relative to modern (meters), set proportional to δ 18 O (Lisiecki and Raymo, 2005) with modern δ 18 O corresponding to 0 m and Last Glacial Maximum δ 18 O corresponding to −125 m.Q 80 is the change in January insolation at 80 • S from modern (W m −2 ) (Laskar et al., 2004).rCO 2 is atmospheric CO 2 in units of preindustrial level (280 ppmv), used mainly for deeper time (pre-Pliocene) experiments.For fixed pre-industrial CO 2 , w c varies between 0 for glacial maxima, 1 for modern-like climates, and 2 for warmest interglacials.w c is converted to 3 weights for those 3 climates (each between 0 and 1, summing to 1): which are used to alter the ocean temperature and basin factor from Eq. ( 18): In order for the model to represent vertical tidewater faces, and to avoid whole grid-cell jumps in the advance and retreat of ice shelves, floating ice is allowed to occupy a subgrid fraction of cell area, f e .This is only applied at ice shelf edges adjacent to open ocean; for interior shelf and all grounded points, f e = 1.The motivation and method here closely follow Albrecht et al. (2011) for the PISM-PIK model.
For floating ice cells adjacent to open ocean, the sub-grid actual thickness (within the area f e ) is estimated based on the thickness of adjacent, presumably upstream, ice (Albrecht et al., 2011).All adjacent points are examined, and the maximum of their ice thicknesses (h) is taken, but only if they are grounded, or are floating and not themselves adjacent to open ocean.Furthermore, if grounded, the interpolated thickness at the grounding line is used.
Then this maximum thickness, h max (m) say, is reduced to allow for typical downstream thinning into the cell in question: h e = max h max max (0.25, e (−h max/ 100) ), 30, h where the minimum of 30 m avoids very thin shelves, and h e can also not be less than the current cell-mean thickness h.h e is the estimated actual ice thickness within areal fraction f e of the cell in question.Finally, to implicitly conserve ice mass, the fractional area occupied by ice in this cell is reset to where h is the cell-mean thickness h (ice volume divided by cell area).Note that the settings above are only done for floating ice cells adjacent to open ocean, otherwise f e = 1 and h e = h.The variable f e is used elsewhere in the model to scale quantities that truly depend on area of ice, i.e., surface mass balance and oceanic melting are both multiplied by f e in the ice thickness evolution as in Eq. ( 14).Also, as mentioned in Sect.2.2, no advective flow of ice is allowed out of a cell with f e < 1.

Calving at ice-shelf edge
There has been considerable recent activity in modeling calving of tidewater glaciers and ice shelves, in part because the extent of floating ice can affect the amount of back stress (buttressing) at the grounding line, and hence the stability of grounded ice upstream (Scambos et al., 2004).Various mechanisms or triggers have been represented in models, including ice thickness over flotation, penetration of crevasses and surface water, and large-scale stress fields (reviewed by Benn et al., 2007; also for instance Alley et al., 2008;Nick et al., 2010;Levermann et al., 2012), but there is little consensus on the main mechanism or mechanisms.
The calving parameterization here is based on the largescale stress field, represented by the horizontal divergence of floating ice velocities.It shares the same motivation as earlier studies by Doake et al. (1998) and is similar to the parameterization in PISM-PIK (Martin et al., 2011;Winkelmann et al., 2011;Levermann et al., 2012), but without using principal strains, i.e., with no distinction between along-flow and across-flow strains, as in Amundson and Truffer (2010).Inclusion of fracture propagation (e.g., Hulbe et al., 2010;Albrecht and Levermann, 2012), multiple stable states (Levermann et al., 2012) and other calving mechanisms are deferred to future work.
First, the divergence of floating ice shelf points div is calculated as using ū and v from the solution of the SSA Eqs.(2a, b) above.This is done only for floating grid points with full fractional cover (f e = 1, Sect.2.9), and propagated by nearest-neighbor value to those on the shelf edge with f e < 1.Then, for points at the shelf edge adjacent to open ocean, the grid-mean calving loss CMB (m a −1 , used in Eq. 14) is set as a weight between two values: where the weight w c = min (1, h e /200).Here, h e is the subgrid thickness of ice within fraction f e (Sect.2.9), and dx is the grid cell size.All units are meters and years.For thin shelves (h e << 200 m), calving is simply weighted towards a constant value of 30 m a −1 .For thicker shelves, it is weighted towards a value proportional to divergence div (a −1 ), but only for positive div.
The thickness h e and grid size dx enter in Eq. ( 26) because 3×10 5 max (div, 0) represents the calving rate (i.e., average horizontal speed of erosion of the shelf edge into the interior, U c in Benn et al., 2007), but CMB here is the rate of volume of ice removed from the cell divided by cell area, so the expression is multiplied by h e dx/dx 2 .
The magnitude of the 3×10 5 coefficient (m) is reasonable on scaling grounds.For a steady-state edge position, the calving rate (U c = 3×10 5 div) must balance the advective ice velocity just upstream of the edge (U T ).For the large West Antarctic shelves, ice velocities change significantly upstream on scales of several 100s km, L T say, so the divergence at the edge is on the order of U T divided by L T .In that case, the parameterized U c = 3×10 5 U T / L T , which is the same order as U T as required for steady state.
CMB is further modified for seaways mostly surrounded by land, represented by the angle subtended to open ocean, A a .This quantity is also used to modify oceanic melt (Sect.2.8, Eq. 18).As discussed in that section, these modifications are needed to allow regrowth of thin shelves in central West Antarctic seaways following a major WAIS collapse (in contrast to the vigorous calving at the edges of the thicker Ross and Weddell shelves today).It can be motivated by considering the effects of icebergs clogging in the restricted seaways, possibly creating a melange that inhibits further calving, but this needs to be explored by future modeling (cf.Vaughan et al., 2011).The calving loss rate CMB is reduced by The divergence div and calving loss given by Eqs. ( 26) and ( 27) are shown in Fig. 5 for a modern nested West Antarctic simulation.In practice, ice-edge thicknesses are often considerably less than 200 m, so the weight w c in Eq. ( 26) is ∼0, and CMB is close to the constant 30 m a −1 in many regions, as seen in Fig. 5.This will be improved in a new calving parameterization under development (see below).For past climates, calving is reduced for cooler environments, similarly to ocean melt in Sect.2.8.This is somewhat ad hoc, because the dependence of divergence on calving does not directly depend on temperature, as some of the other mechanisms mentioned above.But we find that calving must be reduced in order to allow grounding lines to expand as observed during glacial maximum periods.
where w lgm , w mod and w hot are the 3 climate weights corresponding to glacial maxima, modern-like and warm interglacial conditions (as in Sect.2.8, Eq. 20).We are currently developing a new calving parameterization with surface melt dependence, which may avoid the questionable dependencies in Eqs. ( 27) and (28).

Oceanic melt at vertical faces
The parameterization of sub-grid areal fraction in Sect.2.9 allows tall vertical ice faces to be in contact with the ocean, including tidewater fronts extending one grid cell from deep grounding lines.Observations at Greenland calving faces show that oceanic melting of the submerged ice front can be up to a few meters per day (Rignot et al., 2010).A parameterization of the actual circulation and melt rates at a vertical face (Motyka et al., 2003) is not yet in the model.As a placeholder for now, we calculate the area of each vertical face in contact with the ocean, and simply apply oceanic sub-ice melt rates from Sect.2.8 to that area.For any ice cell adjacent to and in contact with open ocean, the vertical extent of submerged ice is z = ρ i ρ w h e for floating ice (29a) where S is sea level and h b is bed elevation.For each of the (up to 4) neighboring cells with no ice and open ocean, z is multiplied by the length of the interface (dx for Cartesian grids) and by that cell's oceanic sub-ice melt (OMB from Sect.2.8).These are summed, and divided by the cell area (dx 2 ) to yield the cell-mean loss of ice due to face melting FMB used in Eq. ( 14).

Bedrock deformation
As in Huybrechts and de Wolde (1999) and Ritz et al. (1997Ritz et al. ( , 2001)), the response of the bedrock to the changing ice and ocean load is a combination of time-lagged asthenospheric relaxation towards isostatic equilibrium, and modification of the applied load by the elastic lithosphere.The treatment here exactly follows Huybrechts and de Wolde (1999).The downward deflection w b of the fully relaxed response (as if the asthenosphere had no lag) is given by where D = 10 25 N m is the flexural rigidity of the lithosphere, ρ b is the bedrock (asthenospheric) density and g is gravitational acceleration.A lower value of D (∼10 23 to 10 24 N m) can optionally be used for West Antarctica (cf.Stern and ten Brink, 1989).The applied load q is q = ρ i gh + ρ w gh w − ρ i gh where h is ice thickness, h w is ocean column thickness, and h eq and h eq w are their values in the equilibrium state (see below).
Equation ( 30) is solved by a Green's function method.The response to a point load P (q times area) versus radial distance r is where kei is a Kelvin function of zeroth order (Brotchie and Silvester, 1969), and L = (D/ρ b g) 1/4 = 132 km is a flexural length scale.w p has significant amplitude within several Llengths of the point load.The w p are summed over the individual point loads of all grid cells (with P = q× cell area) to give w b (x, y), the deflection of the bedrock surface from equilibrium that would occur if the asthenosphere relaxed instantaneously.This is assumed to be proportional to the unbalanced pressure at the top of the asthenosphere due to the load alone (Brotchie and Silvester, 1969).The actual bedrock rate of change is given by where h b is current bedrock elevation, h eq b is its equilibrium value, and τ is 3000 yr.
The equilibrium state (h eq and h eq w in Eq.31, h eq b in Eq. 33) is taken to be the modern observed, assuming that any glacial isostatic adjustments still to occur from the last deglaciation can be neglected (cf.PD12 Appendix B).Equivalently, at the start of a run, the bedrock model alone can be spun up for several 10 000s years with all ice removed, and the resulting ice-free equilibrated state can be used to define h eq b , h eq w (and h eq = 0).

Time steps, adaptive time stepping
The main ice-dynamical time step t i (for Eq. 14) is selected for most experiments depending on model resolution, for instance ∼0.1 to 0.5 yr for 5 to 10 km, ∼0.5 to 1 yr for 20 km, and 2 to 5 yr for 40 km.There is an option for adaptive time stepping that circumvents numerical instabilities, as follows.A restart file is saved at regular time points during a run (spaced ∼1000 yr apart typically).If a numerical exception (NaN) occurs or if physically unreasonable values of ice thickness, temperature or velocity are detected, the simulation reverts to the previous time point using the last restart file, and tries again to run through the next 1000 yr with the time step halved.If an anomaly still occurs during the next 1000 yr, the process is repeated, and is attempted up to 4 times (i.e., with time steps as small as 2 −4 × the nominal value) before aborting.If an attempt makes it through the next 1000 yr successfully, the time step is reset to the nominal value and the run continues on.
For the NetCDF history files, no special action is needed if this adaptive time-looping occurs, because the model snapshots have a unique time index and overwrite any previous snapshot with the same time value.For sequential (ascii) files that contain output at regular intervals, marker records are written that allow a postprocessing program to recognize any time-looping and delete repeated sections as needed.The adaptive-time-stepping capability can be convenient near the start of experiments that are initialized to a state far from equilibrium with the boundary conditions (e.g., modern ice sheet and other geologic time periods).In those cases, blowups and adaptive time looping tend to occur in the first few hundred years, after which the model becomes adjusted to the boundary conditions and the run continues normally.
Other components of the model are time-stepped or reset at greater intervals.The various intervals are as follows: -Ice thickness and dynamics Eq. ( 14): t i , depends on resolution as above.
-Ice and bed temperatures Eq. ( 15): 50 yr, or t i for 10 km resolution or less.
-Recalculating mass balance on ice surface (Sect.3): 50 to 100 yr.At intervening times, recalculation is done for any ice points whose elevation changes by more than 50 m.

Input datasets and climate forcing
Modern Antarctic input fields are obtained mainly from the ALBMAP v1 dataset at 5 km resolution (Le Brocq et al., 2010).The fields used to determine the equilibrium ice and bedrock state discussed in Sect.2.12, with ALBMAP names in parentheses, are ice surface elevation (usrf), bedrock topography (lsrf, topg), and ice thickness (thk).
Various geothermal heat flux maps can be used in the model (Shapiro and Ritzwoller, 2004 (bheatflx shapiro);Fox Maule et al., 2005 (bheatflx fox); Pollard et al., 2005), but these differ considerably from each other on regional scales with noticeable effects on modern results (see next section).Rather than choose one or the other, in the nominal model we specify a simple two-value pattern, with 54.6 mW m −2 under EAIS and 70 mW m −2 under WAIS.
For runs with parameterized climate, observed annual accumulation rate P (van de Berg et al., 2006 (accr)) and surface air temperature T a (Comiso, 2000 (temp)) are used to calculate modern surface mass budgets, as follows: 1. First, T a and P are horizontally interpolated to the ice model grid, and vertical lapse rate corrections are applied: where γ = .0080 h s is the model surface elevation and h obs s is the modern observed elevation interpolated to the ice grid (similarly to Huybrechts, 1998;Ritz et al., 2001).
2. A sinusoidal seasonal cycle is added to T a , giving monthly air temperatures with peak-to-peak amplitude 20 • C at sea level, increasing linearly with elevation to 30 • C at 3000 m, and 30 • C above (based roughly on GCM climates in the GENESIS v3 model).3. A basic positive degree -day (PDD) scheme is applied to the monthly cycle, with coefficient .005m of melt per degree day.Monthly precipitation P is either rain or snow depending on whether monthly air temperature is above or below 0 • C. Any melt or rain immediately becomes mobile and percolates into the ice sheet (Sect.2.7).For modern runs, there is very little surface melt or rain on Antarctic ice.For paleo and future runs with significant melt and rain, a more detailed PDD scheme is available with seasonal refreezing, snow with liquid storage, distinct snow vs. ice PDD coefficients, and allowance for diurnal and synoptic variability (cf.Marshall et al., 2004).In future work we plan to include insolation explicitly (van de Berg et al., 2011).
4. The surface ice temperature, needed as a boundary condition in Sect.2.7, is assumed to be the annual mean of min [monthly air temperature, 0 Instead of parameterizing climate, the model can be driven by a global or regional climate model (GCM or RCM).The climate model provides monthly air temperature and precipitation to the interpolation and PDD schemes in steps 1 and 3 above (e.g., DeConto and Pollard, 2003;Koenig et al., 2011;DeConto et al., 2011), or provides its own annual surface mass budgets calculated with full climate-model physics directly to the ice model.

Modern results
In this section, some basic model results for present-day Antarctica are compared with observations.These simulations have been run to equilibrium with the modern climate, so the comparison ignores any remaining glacial isostatic adjustments in the real world, which are relatively small compared to modern biases (PD12).As discussed in PD12, further work is planned with transient runs through the last deglaciation and extensive comparisons with paleo data (following Briggs et al., 2011;Whitehouse et al., 2012).
Figure 6 compares ice surface elevations h s with observed, using the model with parameterized modern climate (Sect.3) and inverse-derived basal sliding coefficients C(x, y) (Sect.2.4; PD12).Due mainly to the inverse-derived C(x, y), model elevations are within a few 10s meters of observed in most regions.Over the Transantarctics and some other mountain ranges, there are small patches with eleva-tions a few hundred meters too high.As discussed in PD12, these are thought to be due to insufficient sliding through deep troughs cutting through the mountains, only partially compensated by the sub-grid topographic parameterization in Eq. ( 12); however, further work is needed to test that hypothesis.
Much the same level of accuracy in h s is maintained at different resolutions (20 km and 40 km in Fig. 6; 10 km nested in PD12), which is somewhat surprising for regions such as the Siple Coast with ice streams that are scarcely resolved at 40 km.Apparently the proto-streaming at 20 and 40 km does capture basic features such as interleaved unfrozen vs. frozen beds (Fig. 6d, g), and provides the correct overall flux to the grounding line.(At 10 km resolution, individual Siple Coast ice streams are simulated quite realistically, including century time-scale rerouting and stagnating; PD09 Supp.Inf.).
The model grounding line positions, ice shelf thicknesses and extents are the combined result of the model's SSA and SIA dynamics, grounding line flux prescription, and sub-ice oceanic melt, calving and sub-grid pinning parameterizations described above.They are not completely independent of model resolution (Fig. 6), but the effects of resolution are minor and considerably smaller than other model uncertainties.The major Ross and Filchner-Ronne shelves and grounding lines are reasonably realistic, except that the Ronne grounding line has retreated about 200 km too far south (roughly between the Ellsworth Mountains and the Foundation Ice Stream), causing a pronounced low patch in Fig. 6c and f.Other smaller-scale grounding-line errors are seen in Pine Island Bay, Lambert Graben, and especially on the western Peninsula where George VI Sound (between Alexander Island and the mainland) is overridden with thick grounded ice in the model.The latter errors may require higher-resolution modeling and/or coupling with ocean models to correct entirely, but apart from George VI Sound, the errors are not huge and basic regional features are captured.
Figure 7 examines the model sensitivity to the prescribed geothermal heat flux (GHF) map.As noted above, GHF datasets differ significantly even at large scales (Fig. 7, top row).In PD12 we found that these differences can be accommodated by small adjustments in the inverse-derived distribution of basal sliding coefficients C(x, y) (see Sect. 2.4).Here, we show the sensitivity of the model with fixed C(x, y) used in this paper.As shown in Fig. 7, the various GHF maps cause regional and small-scale differences in surface ice elevation of ∼ 100 to 200 m, and significant changes in basal freezing vs. melting patterns.Analogous results are described in Pattyn et al. (2010) for Antarctica, and Rogozhina et al. (2012) for Greenland.
Modern bedrock elevations are also quite close to observed over most regions, showing that the bedrock model in Sect.2.12 is reasonably realistic (Fig. 8).The largest dif-ferences are caused by two main grounding-line errors mentioned above, on the Ronne coast and George VI Sound.However, as discussed in PD12 (Appendices D, E), some of the general agreement may be fortuitous because the model has not taken transient residuals from the last deglaciation into account.
The recent all-Antarctic dataset of surface velocities (Rignot et al., 2011) provides the opportunity to comprehensively test the model velocity field, as shown in Fig. 9 where the dataset (900 m spacing) has been regridded by simple areaaveraging to the model's 20 km grid.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 flow between them.The model's marginal ice thicknesses are generally close to observed (Fig. 6), 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. 9 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.This type of fluctuation could be stochastic in nature (Payne, 1999;PD09 Supp. Inf.).

Past Myr results, and last deglaciation issues
In previous Antarctic applications, the model simulated glacial-interglacial variations in the Pliocene and Pleistocene that are in basic accord with observations (PD09).This includes reasonable first-order agreement with grounding-line retreat during the last deglaciation (last ∼20 kyrs) (PD09 Supp.Inf.), and with surface elevation histories deduced from field data in the Ohio Range (Ackert et al., 2011;Mukhopadhyay et al., 2012).
The applications in PD09 and Ackert et al. (2011) used an earlier model version with a simple two-value specification of basal sliding coefficients (10 −10 and 10 −6 m a −1 Pa −2 , Sect.2.4, Fig. 3a).We repeated the simulations over the last 5 Myrs in PD09 with the current model including the new inverse-derived C(x, y) distribution, and with C(x, y) = 10 −5 m a −1 Pa −2 (maximum value, slippery sediment) specified for all modern ocean beds.As shown in Fig. 10, the results have the same basic features, with predominantly collapsed WAIS in the warm Pliocene, transitioning to larger glacial-interglacial cycles in the Pleistocene, and rarer WAIS collapses in a few Pleistocene interglacials.
However, the maximum glacial ice volumes are less in the current model, ∼29 × 10 6 km 3 compared to ∼33 ×10 6 km 3 in PD09.This is probably due to greater extents of slipperier beds in the new version (PD09 Supp.Inf.Fig. S7A showed much the same effect).Consequently, the equivalent eustatic sea level rise predicted between the Last Glacial Maximum (LGM) and today is only +1.6 m in the new version (3.5 from WAIS, −1.9 from EAIS due to increasing snowfall), compared to +12 m in PD09.These differences in LGM volume are mainly due to thinner ice on the continental shelves and parts of West Antarctica in the current model (Fig. 11).The different ice thicknesses around the margins affect the timing of grounding-line retreat in the major embayments, and simulated relative sea level curves.As mentioned above, we plan to address these issues in upcoming work with transient simulations and model-data comparisons through the last deglaciation.One focus will be the best-fit values of basal sliding coefficients on the continental shelves (cf.Whitehouse et al., 2012).

Conclusions
This paper has described the formulation of a 3-D ice sheet-shelf model, and presented basic validation vs. modern Antarctica.Ice dynamics in the model uses a hybrid combination of the scaled SSA and SIA equations.A parameterization of ice flux across grounding lines (Schoof, 2007) allows grounding-line migration to be captured well, even with coarse (10 to 40 km) grid resolutions.Dynamical tests vs. higher-order models will continue to be important to verify grounding-line behavior, as the model is applied to different domains and scenarios.The model can feasibly be run on continental scales and million-year time scales.Its modern Antarctic ice distributions are reasonably realistic, due in part to an inverse-derived distribution of basal sliding coefficients (PD12).
Although the current parameterizations of sub-ice-shelf melting and calving around Antarctica yield reasonable modern and paleoclimatic results, some aspects are ad hoc and not well constrained by underlying physics.Planned future work includes improvements in sub-ice-shelf oceanic melt, coupling with regional ocean models, and exploration of iceshelf calving parameterizations that depend on surface melt and iceberg clogging that may be needed for past regrowth after marine collapses of the West Antarctic Ice Sheet.
LHS y Left-hand sides of Eqs.2a, 2b (Pa) τ xx Along-flow longitudinal stress at grounding line (Pa) τ f Non-buttressed longitudinal stress at grounding line (Pa) h Ice thickness (m) h s Ice surface elevation (m) h b Bedrock elevation (m) h w Ocean column thickness (m) h eq Ice thickness in bed-equilibrium state (m) h eq b Bedrock elevation in bed-equilibrium state (m) h eq w Ocean column thickness in bed-equilibrium state (m) f e Sub-grid cell-area fraction with ice (0 to1) h e Sub-grid ice thickness within cell-area fraction f e (m) h g Ice thickness at grounding line (m) T Ice temperature ( • C) T m Ice pressure-melting point ( • C) T Homologous ice temperature (relative to pressure-melting point) ( • C) T b w b Lithospheric deflection (m) D Lithospheric flexural rigidity (10 25 N m) L Lithospheric flexural length scale (D/ρ b g) 1/4 (=1.32 × 10 5 m) τ Asthenospheric isostatic relaxation time scale (3000 a) SMB Surface mass balance (m a −1 ) BMB Basal melting of grounded ice (m a −1 ) OMB Sub-ice-shelf oceanic melting (m a −1 ) CMB Calving loss (m a −1 ) FMB Loss due to oceanic melting at vertical faces (m a −1 ) T o Ocean temperature ( • C) T f Ocean freezing point ( • C) K T Transfer coefficient for sub-ice oceanic melting (15.77 m a −1 K −1 ) K Additional O(1) coefficient for sub-ice oceanic melting A a Subtended arc to open ocean (degrees) S Sea level relative to modern (m) T a Annual mean air temperature ( • C)

Fig. 1 .
Fig. 1.Finite-difference staggered grids in the ice sheet-shelf model.h denotes the centers of h-grid boxes, where ice thickness, ice temperatures, and bedrock elevations are calculated.u and v denote the staggered grid points where horizontal velocity components are calculated.
h s = S + h (1−ρ i /ρ w ).(S here is sea level, and h b is bed elevation).

Fig. 2 .
Fig. 2. Idealized flow-line model tests, similar to basic MISMIP (Pattyn et al., 2012), with uniform surface mass balance, an ice divide at the left-hand boundary, a forward-sloping bed into ocean, and using surface-mass-balance increments to q g (see text).(a) Geometry showing sloping bed and ice sheet profiles.(b) Model equilibrated grounding-line positions vs. 1/rheological coefficient A, for various grid sizes and initial states.Solid line shows the analytic solution (Schoof, 2007).(c) As (b) except showing model error (model minus analytic grounding-line position), divided by grid size.

Fig. 4 .
Fig. 4. (a) Sectors used in sub-ice oceanic melt parameterization.Yellow: Amundsen and Bellingshausen Seas, and western Peninsula.Blue: Weddell embayment.Purple: East Antarctica.Red: Ross embayment.(b) Sub-ice oceanic melt rates (m a −1 ) in modern simulation with 20 km resolution.The average values for each major shelf are reasonable(Nicholls et al., 2009;Olbers and Hellmer, 2010;Dinniman et al., 2011), although somewhat lower for the Ross.Rates are noticeably larger nearer the grounding lines due to the depth dependence of the freezing point T f in Eq. (17), especially in Pine Island and Prydz Bays, but not noticeably for the flatter Ross Ice Shelf.
z = 1 is equal to the vertical conductive flux at the top of the bedrock (see below) plus any basal shear heating Q b = τb .ũb ; or, if T(x,y,1,t) would exceed the basal pressure melt point T m , then it is set equal to T m and the imbalance in conductive fluxes plus Q b is used to melt basal ice.For floating ice, the basal boundary condition is simply T(x,y,1,t) = T m .(Oceanic melt rates under floating ice are parameterized separately in Sect.2.8).

[
longitude, latitude] = [160 to 180, all] or [−180 to −140, all] or [−140 to −120, < −77].T o = −1.5 K = 1 At this point, T o and K represent conditions under modern exposed shelves.For the West Antarctic sectors, ocean melt is further reduced based on subtended arc to open ocean A a (degrees), i.e., the angle formed by the set of all straight lines from the point in question that reach open ocean without hitting land (as in PD09).
1.7 w lg m + T o w mod + 5 w hot (21a) K = K w lg m + K w mod + 8 w hot (21b) Finally, T o and K are modified for distal locations, to prevent ice shelves from expanding into the Southern oceans far from Antarctica.This is based on ocean bathymetry (h w = sea levelh b ), assuming much warmer waters at depths >∼ 2000 m, with an additional constraint based on the arcto-open-ocean angle A a to ensure this is not done for deep proximal troughs.The final T o and K are used in Eq. (17) in place of T o and K. T o = T o (1 − w d w e ) + max [T o , T dist ] w d w e (22a) K = K (1 − w d w e ) + 10 w d w e (22b) where T dist = −0.5wlgm + 5w mod + 8w hot (22c) w d = max[0, min[1, (h w − 1900)/200]] (22d) w e = max[0, min[1, (A a − 150)/20]]

Fig. 9 .
Fig. 9. (a) Observed surface ice velocity (Rignot et al., 2011), averaged here to 20 km model cells, m a −1 .(b) Model surface ice velocity, m a −1 , in modern simulation at 20 km resolution.(c) Model minus observed log 10 (velocity, 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. 10 .
Fig. 10.Time series of total Antarctic ice volume (10 6 km 3 ) over the last 5 million years, in simulations with parameterized climatic and oceanic forcing dependent mainly on deep-sea core δ 18 O, and slightly on austral summer insolation, with 40 km model resolution.(a) Current model, with inverse-derived basal sliding coefficients C(x, y), and value on continental shelves = 10 −5 m a −1 Pa −2 .(b) Earlier model version as in PD09 (their Fig. 3a) with simple twovalued C(x, y) and continental-shelf value = 10 −6 m a −1 Pa −2 .

Table 1 .
Model symbols and nominal values.
18O (as for Eq.19).rCO 2 is atmospheric CO 2 in units of preindustrial level (280 ppmv), assumed to produce a 10 • C warming in the Antarctic region for each CO 2 doubling.Q 80 (W m −2 ) is the change in January insolation from modern at 80 • S. T a is applied on the right-hand side of Eq. (34a) and so also affects precipitation P in Eq. (34b).
• C]. a is applied to air temperatures, mainly determined by deep-sea core δ 18 O and CO 2 , with a minor effect of austral summer insolation (similarly to past ocean melt variations in Sect.2.8, Eq. 19):T a = 10 S 125 + 10 log(rCO 2 ) log(2) + 0.1 Q 80(35)where S (meters) is eustatic sea level relative to modern, proportional to δ