Skip to main content Accessibility help


  • Access
  • Cited by 35



      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Ice-shelf basal channels in a coupled ice/ocean model
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Ice-shelf basal channels in a coupled ice/ocean model
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Ice-shelf basal channels in a coupled ice/ocean model
        Available formats
Export citation


A numerical model for an interacting ice shelf and ocean is presented in which the ice- shelf base exhibits a channelized morphology similar to that observed beneath Petermann Gletscher’s (Greenland) floating ice shelf. Channels are initiated by irregularities in the ice along the grounding line and then enlarged by ocean melting. To a first approximation, spatially variable basal melting seaward of the grounding line acts as a steel-rule die or a stencil, imparting a channelized form to the ice base as it passes by. Ocean circulation in the region of high melt is inertial in the along-channel direction and geostrophically balanced in the transverse direction. Melt rates depend on the wavelength of imposed variations in ice thickness where it enters the shelf, with shorter wavelengths reducing overall melting. Petermann Gletscher’s narrow basal channels may therefore act to preserve the ice shelf against excessive melting. Overall melting in the model increases for a warming of the subsurface water. The same sensitivity holds for very slight cooling, but for cooling of a few tenths of a degree a reorganization of the spatial pattern of melting leads, surprisingly, to catastrophic thinning of the ice shelf 12 km from the grounding line. Subglacial discharge of fresh water along the grounding line increases overall melting. The eventual steady state depends on when discharge is initiated in the transient history of the ice, showing that multiple steady states of the coupled system exist in general.


Observations of several Greenland and Antarctic ice shelves reveal deep inverted channels oriented in the direction of ice flow and running along the underside of the floating ice. For instance, Motyka and others (2011) showed that, before its massive ice loss beginning in the late 1990s, Jakobshavn Isbr^’s (Greenland) floating ice shelf experienced high basal melt rates which likely led to the formation of a single prominent basal channel several hundred meters wide and tens of meters deep. Pine Island Glacier, Antarctica, also terminates with an ice shelf in contact with warm ocean water (Jenkins and others, 2010), and basal channels up to 200m deep are observed there (Payne and others, 2007; Bindschadler and others, 2011), which have a role in steering the flow of water out from underneath the ice shelf (Mankoff and others, 2012). Fricker and others (2009) also found basal channels several hundred meters deep and ~1.5km wide under the Amery Ice Shelf, East Antarctica. Although the channels at the Amery Ice Shelf likely originate at the juncture of neighboring ice streams, as opposed to having a clearly oceanographic origin, the authors suggested these channels may influence the ocean circulation beneath the shelf. Perhaps the clearest examples of pronounced ice-shelf basal channels resulting from ocean interaction are found at Petermann Gletscher, northwest Greenland (80.5° N, 60° W). Observations of highly spatially variable melt rates by Stewart and others (2004) suggest the existence of deep basal channels, and a series of ice thickness measurements by Rignot and Steffen (2008) confirmed the existence of pronounced longitudinal channels in the ice-shelf base which, within 10 km of the grounding line, increase in amplitude up to hundreds of meters in depth.

Increased rates of delivery of ocean heat to Jakobshavn Isbræ in the late 1990s led to the loss of its ice shelf (Holland and others, 2008; Motyka and others, 2011; Hansen and others, 2012) with a consequent acceleration in grounded ice loss (Thomas and others, 2009). Present-day mass loss in Antarctica primarily occurs in regions where ocean- induced thinning of ice shelves lead to reduced buttressing of grounded ice (Pritchard and others, 2012), and it is likely that future contributions to eustaticsea level from the marine- bedded West Antarctic ice sheet will be driven by ocean- induced melting of ice shelves, with atmospheric forcing causing the shifts in ocean circulation which will deliver ocean heat to the ice shelves (Joughin and Alley, 2011). To project future global mean sea-level changes, it is therefore important to understand in detail the coupled ice and ocean dynamics which determine the overall rates and spatial patterns of ice-shelf basal melting.

Our primary aim was to investigate the possibility that buoyancy-driven ocean currents under the Petermann Ice Shelf bring ocean heat into contact with the ice and generate spatial patterns of melting which, in combination with the ice flow, lead to the formation of basal channels. The model results below show that basal channels at the Petermann Ice Shelf are indeed the product of feedbacks relating ice- shelf basal slopes and ocean melting. However, given the presence of warm water, we will show that an ice shelf which develops deep and narrow basal channels near the grounding line experiences less melting than it would if such channels failed to develop and melting occurred on spatially broader scales. In this sense, we find that Petermann’s narrow channels are a concession to warm water which allows the ice shelf to remain viable. Our second major finding is that, in the regime of ocean temperatures and ice inflow velocities we investigated, melt rates are generally so large that ice deformation is ineffective in balancing the thinning tendency of melting. Instead, the thinning tendency of melting is primarily offset by the down-glacier advection of ice by the nearly spatially uniform background velocity of the ice shelf. Other questions we address in the following are the sensitivity of the coupled ice/ocean system to changes in ambient water temperature, changes in rates of subglacial freshwater discharge, changes in glacier speed and various changes to model parameterizations related to ocean turbulence and mixing.

Petermann Gletscher

Although it is not observed to be unsteady to any significant extent (Thomas and others, 2009), we chose to investigate Petermann Gletscher with the expectation that a better understanding of ocean-induced basal melting in this case will shed light on ice-shelf basal channels more generally. In the rest of the paper we therefore restrict our discussion to this 60 km long ice shelf, which is the ocean outlet for ~6% of the Greenland ice sheet (Rignot and Kanagaratnam, 2006).

In Figure 1 a surface elevation map of the Petermann Ice Shelf derived from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery confirms that the basal undulations measured by Rignot and Steffen (2008) constitute a system of channels running along the ice base, which intensify just seaward of the grounding line and then decay in amplitude at distances >20 km from the grounding line. If we refer to the ice which protrudes to the greatest depths as ‘keels’ and the base of the thinnest ice as ‘channel crests’ then we may observe that the typical keel-to-crest ‘channel amplitude’ is ~300m within 20 km of the grounding line. These deep longitudinal channels are incised into ice which would otherwise be at most 500 m thick. Across-shelf sections of ice basal elevation presented by Rignot and Steffen (2008) show that channel amplitudes well exceed the scale of any bedrock topographic irregularities just upstream of the grounding line and the characteristic width of the channels of ~3 km does not correspond directly to any obvious features in the grounded ice. Ocean observations both under the ice shelf (Rignot and Steffen, 2008) and at the mouth of the fjord seaward of the ice shelf’s terminus (Johnson and others, 2011) show that, below ~200m depth, the ice shelf is immersed in relatively warm Atlantic waters which reach Petermann via the Arctic Ocean.

Fig. 1. Relative surface elevation of the Petermann Ice Shelf derived from stereoscopic ASTER imagery, which gives an impression of the basal topography since the ice hs been floating (personal communication from T. Millgate, 2011).

Ice/Ocean Model

In order to investigate basal channel formation we have created a numerical model with which we compute the evolution of an idealized but Petermann-like floating ice mass in contact with a buoyancy-driven ocean mixed layer (an extended plume model) at the top of the water column (Fig. 2). The model was obtained by coupling the Glimmer Community Ice-Sheet Model (Glimmer-CISM) (Rutt and others, 2009; Lemieux and others, 2011; Price and others, 2011) to a buoyancy-driven ocean mixed-layer model (Holland and Feltham, 2006; Holland and others, 2007, 2009).

Fig. 2. Petermann-like model geometry.

Relation to prior models

At present, several ocean models have been modified to include appropriate ocean boundary conditions and thermodynamic interactions between ice and ocean in sub- ice-shelf environments. For instance, Jenkins and Holland (2002) used an isopycnal coordinate model to investigate the buoyancy-driven ocean circulation under the Filchner- Ronne Ice Shelf, and Holland and others (2003) considered the regional significance of the circulation under the Ross Ice Shelf. Thoma and others (2008) used the same ocean model to investigate the variability of warm circumpolar deep water intrusions into Pine Island Bay in the Amundsen Sea. Mueller and others (2012) used a σ-coordinate model to investigate the effect of tides on melting beneath the Larsen C ice shelf. Heimbach and Losch (2012) used the automatically generated adjoint of the z-coordinate MITgcm model to determine the sensitivities of melt rates under the Pine Island Ice Shelf to ocean temperature. All of these models shed new light on ocean circulation in sub-ice-shelf environments, but none attempted to investigate the feedback arising from the mutual influence of ocean circulation, melt rates, ice dynamics and ice geometry, since they did not include dynamic ice shelves. Initial progress in coupling dynamic ice models with dynamic ocean models was made by Grosfeld and Sandhager (2004) who coupled three-dimensional (3-D) ice and ocean models on a 50year coupling cycle, enabling millennial-scale simulations of an idealized Filchner-Ronne Ice Shelf. Our choice to use a simplified ocean model in our ice and ocean system made it computationally feasible to explore many model configurations and to integrate the system to steady state in various configurations while investigating the coupled interaction of ice and ocean.


The Glimmer-CISM ice-sheet model is based on a set of equations first put forward by Blatter (1995) and Pattyn (2003), which we present here for completeness. Let σij denote the deviatoric stress tensor in the ice, the strain-rate tensor, p the pressure, the velocity field, p1 the ice density and gravitational acceleration.

Stress-balance equation

Assuming the ice is in a stress equilibrium, the flow is governed by the Stokes equations


Stresses, σij , are related to strain rates, according to Glen’s constitutive law, which can be written as


with a strain-thinning viscosity given by where is the effective strain rate. The small value, , is added to prevent the viscosity, η, from becoming singular when . In the following we take n = 3 and the Arrhenius coefficient, A, corresponding (Payne and Dongelmans, 1997) to a uniform ice temperature, T1 , which is held constant. In the control run we set T1 = −10°C.

The Higher-Order (Blatter-Pattyn) Model (HOM) refers specifically to the simplification obtained by neglecting the so-called bridging stresses, σ zx and σzy , in the vertical stress balance in Eqn (1) and neglecting the corresponding shearing terms, ∂w/∂x and ∂w/∂y, in the strain-rate tensor in Eqn (2). Detailed validations of these simplifications in the case of a small aspect ratio of the ice are given by Dukowicz and others (2010) and Schoof and Hindmarsh (2010).

This simplification decouples the computation of pressure and vertical velocity from the computation of the horizontal velocity, (u, v). The vertical stress balance yields


and using σ xx + σyy + σzz = 0 we obtain equations



governing the horizontal velocities u and v.

The numerical solution of Eqns (4) and (5) in Glimmer-CISM involves a Picard iteration, which is based on the observation that if η(x, y, z) is fixed, the resulting equations for u and v are linear and can be solved by well-developed iterative methods. Then η can be recalculated and the process iterated until convergence is obtained. We found that imposing a lower limit, η ≥ η0 = 1 0 x 10-3 m2 s-1, allows the Picard iteration to converge in many cases of interest where convergence would otherwise fail. Note that a minimum viscosity corresponds to a maximum effective strain rate, so this numerical strategy is independent of and additional to the regularization involving

Ice thickness equation

The ice is treated as incompressible with p1 = 910 kg m-3 and therefore the ice thickness, H, evolves according to


where is the depth-averaged horizontal velocity, ȧ is the surface accumulation or ablation rate, ṁ is the basal melt rate and K is artificial thickness diffusivity. The appearance of a factor y here reflects our use of a simple homotopy method to throttle melt rates during model spin-up and in certain perturbation experiments. The final term is unconventional and is introduced out of necessity to artificially smooth the ice in the across-shelf direction. A detailed description of these two modifications appears in Appendix A. In the numerical solution of Eqn (6), the advective part uses incremental remapping, a scheme which was applied to sea- ice advection by Lipscomb and Hunke (2004) and made available for use in Glimmer-CISM by the same authors.

Note that, in accordance with its neglect of bridging stresses, Glimmer-CISM automatically assumes ice columns located above sufficiently deep bedrock will float in local hydrostatic equilibrium. Consequently Glimmer-CISM is not suitable for solving the fully general contact problem determining the position of grounding lines. In order to isolate ice-shelf/ocean interaction as a possible mechanism for basal channel formation, and to avoid model limitations associated with grounding-line migration, we have therefore chosen the bedrock topography to be everywhere deep enough to ensure all ice in the domain is floating at all times. Assuming the ice rapidly achieves local hydrostatic equilibrium everywhere throughout the domain, the lower surface of the ice is given by b = −H10), with ρ0 = 1028 kgm-3, a reference ocean density.

Geometry and boundary conditions

Our primary Petermann-like model domain is a 20 km x 40km rectangle completely filled with floating ice. We model only the first 40 km of the Petermann Ice Shelf, since channels appear to be largely eroded downstream of this position, and the ocean circulation is not likely to be correctly described by our buoyancy-driven model farther downstream, where external water masses likely penetrate under the ice (Johnson and others, 2011). Both ice and ocean models are discretized on a regular Cartesian grid with 250 m x 250 m cells.

On the southern boundary, velocity is fixed as (u, v) = (0, v 0) with v 0 = 1kma 1 in the following calculations, unless stated otherwise. Vertical velocity, w, is zero along the bottom of the inflow boundary and otherwise determined by the hydrostatic relation, Eqn (3), along with the condition of flotation. Along the eastern and western lateral walls, there is no normal flow and we impose a lateral shear ∂v/∂x = (1/η)τ0, where τ0 is a fixed yield stress of 25 kPa along the western wall and –25 kPa along the eastern wall. In other words, we assume deformation along the lateral walls is perfectly plastic and we simply choose a reasonable yield stress (Cuffey and Paterson, 2010) in the absence of relevant measurements. The upper surface is stress-free, while the basal surface has no shear stress and normal stress implicitly equal to that required to support the weight of the local ice column in hydrostatic equilibrium.

Finally, at the northern outflow edge the shear stress vanishes and the normal strain rate, following MacAyeal and others (1996), is determined from the balance


between the depth-integrated normal stress in the ice (left side) and depth-integrated water pressure at the ice front (right side). The heights s and b are of the ice surface and base, respectively, as in Figure 2. Note that Eqn (7) is exact, since it is just the statement of the continuity of stress across the ice front. The Blatter-Pattyn equations imply that the absolute normal stress, azz – p, matches the ice overburden weight, yielding

Then Eqn (7) implies

Assuming u and v are nearly depth-independent, which later numerical experiments verify, we obtain the final relation

which is imposed along the northern edge of the domain.

The initial ice thickness is prescribed everywhere and the inflowing ice along the southern edge is given a prescribed thickness, H0(x), which is fixed for each model run. Because no thickness information enters the domain along the lateral walls (due to no normal flow) or along the northern boundary (due to strictly outgoing velocities) no thickness boundary conditions are prescribed there. Following Rignot and Steffen’s (2008) automated weather station observations, we impose a uniform 1.0ma_1 ablation over the entire upper surface of the ice. Basal melt rates, ṁ, are determined by the ocean model.


The ocean model describes a buoyancy-driven mixed layer on top of a deep ambient ocean layer. Specifically, the ambient ocean is taken to be horizontally homogeneous, with temperature and salinity given by fixed vertical profiles, Ta(z) and Sa(z). The ambient ocean is taken to be at rest. The ocean model provides prognostic equations for the mixed- layer thickness, D, depth-averaged salinity, S, temperature, T, and horizontal velocity,

Ocean layer thickness equation

Integrating the incompressibility condition from z = a to z = b yields the ocean mixed-layer thickness equation


where ṁ ė and ḋ represent melting at the ice/ocean interface, and entrainment and detrainment at the ambient ocean interface, respectively. If the ocean layer is gaining mass by entrainment then ė > 0 and ḋ = 0. Otherwise, when the layer is losing mass by detrainment, ḋ > 0and ė = 0. Entrainment and detrainment are represented separately for later convenience. Subglacial discharge of fresh water, s, is treated as a mass source rather than a boundary condition and therefore appears as a volumetric flux, s, which is always zero except possibly near y =0.

Momentum equations

By depth-integrating the rotating and Reynolds-averaged Boussineq equations we obtain momentum equations




where z = a(x, y, t) = b(x, y, t) – D (x, y, t) is the elevation of the interface between the mixed layer and the underlying ambient ocean, f is the Coriolis parameter, Ah is a horizontal eddy viscosity, g’ = ga −ρ)/ρ0 is the reduced gravity determined by the density, ρ, of the mixed layer and the density, ρa, of the ambient water immediately below the mixed layer, and cd =2.5 x 10-3 is a fixed non-dimensional drag coefficient. The terms -ḋU and -ḋV represent the loss of fluid transport from the mixed layer due to detrainment, ḋ. A corresponding term for entrainment, ė, does not appear since ambient water absorbed into the mixed layer by entrainment is assumed to be motionless.

Terms involving g are the result of integrating the horizontal pressure gradient over the thickness of the mixed layer. The pressure gradient in the mixed layer was obtained by assuming the horizontal pressure gradient vanishes in the deep ambient layer and then applying the assumption of hydrostasy to integrate upwards.


The entrainment rate, ė, is given by rearranging the formula


obtained by Gaspar (1988) by considering the turbulent kinetic energy (TKE) budget in each vertical column of the mixed layer. Terms on the left-hand side are sinks of TKE, while production of TKE from all sources minus molecular dissipation of TKE is parameterized by the term on the right-hand side. In this formula, ΔBa ≥ 0 is the buoyancy difference across the mixed layer/ambient ocean interface and ΔBb ≥ 0 is the buoyancy difference between the mixed- layer water and water in contact with the ice base z = b, with temperature and salinity denoted Tb and Sb, respectively The non-dimensional coefficient, μ, is a tuneable parameter which we take to be 2.5 (which is considerably larger than the upper estimate μ = 0.5 suggested by Gaspar, 1988)) and is the friction velocity.

If Eqn (11) predicts a negative value for ė we instead calculate the increment, ΔD, by which the mixed-layer thickness must decrease in order to reduce to the Monin- Obukhov thickness, (2 μU* 3) /(ΔBbṁ). The detrainment rate is then ḋ = ΔD/Δt, where Δt is the ocean time-step. In other words, detrainment is modelled as an instantaneous retreat to the Monin-Obukhov depth in the case where TKE production is inadequate to mix down meltwater produced at the ice/ocean interface. Former mixed-layer water left behind by the retreating ambient ocean interface is treated as though it instantaneously mixes into the ambient water without changing the ambient water’s characteristics.

Scalar quantities

Ocean mixed-layer temperature and salinity are governed by depth-integrated advection-diffusion equations



where T a and S a are the ambient temperature and salinity just below the ambient interface z = a, Kh, is an eddy diffusivity, and γT is an exchange velocity for turbulent heat transfer from the mixed layer into the ice/ocean boundary layer. We always take Kh = Ah, and for the control case below we use Ah = Kh = 10m2s-1, which gives a gridscale Reynolds number which is O(1) nearly everywhere. Subglacial discharge, ṡ, adds fresh water at 0°C. Since this does not change the depth- integrated heat content, DT, or salt content, DS, of the mixed layer, s does not appear directly in Eqn (12) or (13).

Thermodynamic relations

Densities ρ and ρa are determined using a linearized equation of state,


where the coefficients βT = 3.87 x 10 C-1 and βS = 7.86 x 10-4psu-1 are appropriate for T near T0 = -2.0°C and S near S0 = 34.5 psu. Water at the unresolved molecular interface between ice and ocean is assumed to have properties Tb and Sb that are constrained to lie on a linearized freezing temperature line


where α = −5.73 x 10-2°C psu-1, β = 8.32 x 10-2°C and λ = 7.61 x 10-4°Cm-1. Note that the ice draft, b, is negative so T b decreases with depth. Turbulent heat and salt fluxes into the ice/ocean boundary layer are prescribed from bulk formulae and are balanced by melting or refreezing according to (Holland and Jenkins, 2001)



Here ṁp1c1(Tb – T1) represents the energy needed to heat the ice from its internal temperature, T1 , to the in situ melting point, Tb, £ = 3.35 x 105 Jkg-1 is the latent heat of ice fusion, c o = 3984 J kg 1 C 1 is the specific heat capacity of sea water and c l = 2009 J kg 1 C 1 is the specific heat capacity of ice.

The turbulent exchange velocities, YT and YS, are given by



where v0 = 1.95 x 10 ~6 m2 s-1, Pr = 13.8 and Sc = 2432 are, respectively, the molecular viscosity, Prandtl number and Schmidt numbers for sea water.

Ocean boundary conditions and ambient ocean

Free-slip and zero-scalar-flux boundary conditions are prescribed along the west, south and east boundaries of the ocean mixed layer. Along the northern edge we impose ∂U/∂y = ∂V/∂y = ∂S/∂y = ∂T/∂y = 0 as an approximation to an open (radiation) boundary condition with rapid outflow (Durran, 1999). Numerical instability can arise if V is allowed to take on negative values along the northern boundary, so this is not permitted.

Some ocean features inferred from observations in the vicinity of Petermann Gletscher cannot be captured by the model. For instance, cold and fresh fjord surface water may intrude under the ice (Johnson and others, 2011) and the resulting stratification may make it possible for a buoyancy-driven flow to overshoot its neutral density depth and subsequently separate from the ice base (e.g. Baines, 2008). We cannot model this separation with the equations presented above, so we always select Ta(z) and Sa(s) in a way that makes the ambient water neutrally stable. Since the ocean mixed layer is always a mixture of ambient water and fresh meltwater, it is lighter than ambient water at any depth and it is therefore impossible for it to overshoot its neutral density level.

Remark on forcing of the ocean

It is easily verified that the ocean equations admit a degenerate solution, D = U = V = ė = = 0 and (T, S) = (Tb, Sb). Numerical experimentation shows that without appropriate forcing on the right-hand side of the thickness equation, Eqn (8), the mixed layer will tend toward this degenerate solution through a combination of detrainment and velocity divergence. In particular, note that if speed, is very small, melting causes the layer to detrain, according to Eqn (11), to a correspondingly small thickness. Once this occurs, terms on the right-hand sides of Eqns (9) and (10) are also small, unless ∇ρ or ∇a are large. Ultimately, the ocean layer needs to be continuously fed from upstream in order to persist. This requires either a boundary condition on the plume thickness (as in Holland and Feltham, 2006; Holland and others, 2007, 2009) or some type of forcing to prevent the layer tending to zero thickness.

We experimented with three different methods for driving the mixed-layer circulation. The first approach is simply to artificially increase the entrainment rate, e, as necessary to maintain a minimum layer thickness of D 0 everywhere in the ocean domain. Physically, this is like postulating a source of TKE not otherwise accounted for in Eqn (11), such as shear due to tidal oscillations or flow in the ambient layer (which we do not model), that enhances entrainment when the mixed-layer depth is small. With the choice D0 = 1.0 m, used in all the following simulations unless otherwise indicated, the regions where artificial entrainment is needed to maintain this minimum thickness are typically along the deep ice keels and along the boundaries, places where mixed-layer velocities tend to be divergent.

A second approach, again invoking tides, is to replace the formula with the formula


where U2 *0 represents a constant drag on the underside of the ice due to unresolved (i.e. tidal) flow. This increases the amount of TKE available for entrainment and deepens the mixed layer, driving circulation. The form of Eqn (20) arises from supposing that the sub-shelf circulation, due to tides is uncorrelated with the strictly buoyancy-driven circulation, over a tidal cycle, so that everywhere, with the average being taken over a tidal cycle at a fixed location. Then the friction velocity, U*, on the timescale of a tidal cycle or longer would be better estimated as

where (|U0|) represents the rms speed of the unresolved tidal flow. For all the calculations presented below, unless otherwise indicated, we use U*0 = 0.0025 ms-1, corresponding to an rms tidal speed of 0.05 ms-1, given our choice of Cd. Note that increasing U* in this way tends to increase entrainment rates as well as turbulent heat and salt transfer into the ice/ocean interface.

The third approach neglects tides and supposes the circulation is maintained by the discharge of buoyant fresh water along the grounding line. Specifically, in this approach certain patterns of subglacial discharge, s(x, y), are imposed along the southern edge of the domain where y < 1.0 km. This strongly increases ∂p/∂y and ∂a/∂y near the inflow edge, causing the mixed layer to accelerate and entrain ambient water. We use s = 0 in the control case below, but experiment with positive s in a set of perturbation experiments described below.


Both models use the same uniform rectangular grid with a grid spacing of 250 m for all calculations presented below.

The two models exchange information via the melt rate, m, and the ice draft, b. For example, positive melt rates m tend to thin the ice and add mass and buoyancy to the mixed layer. In turn, the ice draft, b, influences the ocean momentum equations through the ‘interfacial-flattening’ terms involving ∇a = ∇b —D.

In general, for ice geometries encountered in our experiments, if the ocean model is allowed to run indefinitely against a fixed ice geometry it reaches an essentially steady state within one model day, which is consistent with flushing of the 40 km domain by model currents with a velocity scale of ~5 cms-1. In certain regions, however, especially where the mixed layer is thin and fast with strong gradients in thickness, oscillations in the thickness and velocity with a period of ~3 hours prevent melt rates from becoming steady in those regions. Whenever these oscillations arise they always propagate at an angle of 45° and have a checkered spatial pattern, which implies they are spurious numerical errors in nature.

For all of the following coupled experiments, the ocean is initially run to a steady state with respect to the initial ice geometry. During each ice time-step the ocean is allowed to run for a subcycling time of 0.1 days, using a 60.0 s time-step. Averaging the melt rates over this subcycling time window partially filters out the spurious oscillations mentioned above and this melt rate is passed to the ice model to be applied over the ice time-step.

Since the ice time-step is At = 0.05 years, the ocean adjusts for a total of one model day for each half-year of ice evolution. Changes in the ice fields are generally small over this timescale and therefore the ocean is never very far from a steady state. This coupling approach permits the model ocean to be integrated over the long timescales associated with the ice dynamics.

Model Results

Petermann-like control case

Here we describe a model run that reproduces several salient characteristics of Petermann Gletscher, including deep basal channels.

Choice of inflow ice profile, H0(x)

A thickness profile H0(x) must be imposed as a boundary condition along the southern edge of the domain. We experimented with many possibilities before selecting

for the control case. We refer to this profile, which has a mean thickness of 600 m with four primary 50 m undulations and twelve secondary 25 m undulations superimposed on it, as the ‘4 plus 12’ profile (Fig. 3b).

Fig. 3. (a) Fixed melt rate applied to the ice to achieve the pre-control steady state. (b) Ice draft for this steady state.

While this choice is somewhat arbitrary, we justify our use of it by the fact that simpler choices failed to develop Petermann-like basal channels or led to excessive basal melting and a consequent failure to achieve a coupled steady state. For example, we found that a flat profile, Ho(x) = 600.0, led to excessive melting along the eastern boundary of the shelf (melting completely through the ice) as the ocean layer formed a jet banked against the eastern boundary under the influence of the Coriolis force. As we show below (Fig. 12; Table 1), profiles of the form lead to an excess of melting. For k > 7 the model fails to produce the strikingly deep (200–300 m) channels found under the Petermann Ice Shelf. Using k = 6 or 7 would be a reasonable control case, but the resulting steady states bear somewhat less resemblance to observations than the ‘4 plus 12’ choice, given that three or four major channels appear to emerge in the Petermann Ice Shelf once the shorter-wavelength variations die away (Fig. 1). CReSIS (Center for Remote Sensing of the Ice Sheets, University of Kansas) bed elevation data compiled from 1995 to 2011 show that, upstream of the grounding line, undulations in the bed shape ~50 m in depth with wavelengths of 2–6km are typical (data downloaded from in November 2011). We therefore proceed with the ‘4 plus 12’ thickness profile, as a case which exhibits a degree of realism.

Table 1. Time (years) at end of simulation for non-steady cases in Figure 12

Spin-up procedure

To spin-up the control case we imposed the inflow ice thickness profile, H0(x), described above, and ran the ice component to a steady state with an imposed basal melt field uniform in the x-direction and given by a piecewise linear function in the y-direction approximating melt rates inferred from observations by Rignot and Steffen (2008). The resulting steady-state ice draft and the imposed basal melt rates are shown in Figure 3. We refer to this ice geometry as the ‘pre-control’ case. The absence of deep basal channels in this steady state allows us to rule out a purely glaciological origin for the observed basal channels. On the contrary, ice deformation tends to very slightly smooth out the cross-shelf irregularities as the ice moves north, if only by a few tens of centimeters per year (i.e. not enough to remove the incipient channels already present in the inflow profile). The artificial thickness diffusivity, K, was set to zero for this phase of the spin-up.

Next, coupling of the ice to the ocean model was activated and the fully coupled model was allowed to run from the pre-control state until a coupled steady state was obtained. The ambient ocean temperature and salinity profiles used in the control case are shown in Figure 4. The temperature profile is a piecewise linear function approximating the ocean temperatures measured below the Petermann Ice Shelf (Rignot and Steffen, 2008). The salinity profile was chosen to match the observed salinity at depth in the ice-shelf cavity (34.75 psu) and, to prevent mixed-layer separation as explained above, has vertical variations tied to the temperature variations so as to yield a water column of uniform density.

Fig. 4. Fixed temperature and salinity profiles used to represent ambient water in the ocean model. The blue curve is temperature and the red curve is salinity. Density is uniform throughout the ambient water column.

Overview of the control case results

Figures 5 and 6 give an overview of the model results and show the essential relationship between ice and ocean that maintains deep basal channels in the steady state. In Figure 5, velocity vectors for the ocean layer show that water circulates especially vigorously near the southern edge, where it banks against and flows up the steep channel flanks and headwalls. The areas where the ocean layer flows the fastest are precisely the areas where the melt rates are the largest (also see Fig. 9). The spatially variable melting across the shelf removes ice from the base of the shelf as it passes by, cutting out a profile that is then extruded into channels by ice advection. In Figure 5c we show the ‘channel depth’, which we define, for ice in a given crossshelf section, to be the height of the ice-shelf base relative to the deepest ice along that section. Regions with channel depths of up to 250 m are found just north of the high- intensity melt regions. Figure 6 shows sections of the ice and ocean layer in the region where channels develop. Channels with wavenumber k = 4 are amplified significantly with the shape of the channel profile becoming sharper between y = 1 and y = 9 km. Channel deepening is biased by the Coriolis-induced banking of water against the eastern flanks of basal concavities. This results in the elimination of k =12 concavities at y = 0.8, 5.8, 10.8 and 15.8 km and the amplification of the k = 4 channels and the k =12 concavities at y = 4.2, 9.2, 14.2 and 19.2 km. In the k = 4 channels, melting is slightly stronger to the east, which makes the eastern channel flanks steeper than the western flanks and causes the channel crests to migrate slightly to the east. By y = 40 km the channel amplitudes have decreased significantly and the profile has become much smoother. The intensification of basal channels downstream of the grounding line and the moderation of channels further along the ice shelf are notably similar to the character of the ice revealed by both radar observations and satellite altimetry (Rignot and Steffen, 2008; personal communication from T. Millgate, 2011<). Modeled ice velocities (Fig. 7) are found to be much less spatially variable and indeed deviate relatively little from the ‘background velocity’, (0, v0), imposed as a boundary condition on ice entering across the southern edge.

Fig. 5. (a) Ice draft variations are due to intense melting near the southern inflow edge, where the ocean mixed-layer circulation (velocity vectors shown) is vigorous. (b) Melt rates are strongly varying in the x-direction. (c) Draft anomalies relative to the deepest ice at each y location show the presence of channels up to 250 m deep.

Fig. 6. Cross-shelf sections showing the ice base, z = b(x, y), in blue and the bottom of the mixed layer, z = a(x, y), in red. The gap between the curves in each pair shows the thickness of the mixed layer. Note the accumulation of ocean mixed-layer water in the crests of channels, with deeper water to the right of the crests, due to deflection by the Coriolis force.

Fig. 7. Colors show the draft of the ice base. The plotted arrows show (u, v − V0), the horizontal ice velocities with the inflow velocity V0 = 1 km a-1, removed. These velocities, which are associated with ice deformation, are generally <100 m a-1 and lack significant spatial variation on the width scale of basal channels.

Fig. 8. At each location, y, along the ice shelf, the Fourier transform, Ĥ(k,y), of the ice thickness variations, , was computed. Here we plot log where k, the wavenumber of variations in the x-direction, is plotted along the horizontal axis and y is on the vertical axis. This ‘cross-shelf spectral density’ plot shows that variations in H0(x) along the southern edge of the domain with wave numbers k = 4 and k =12 are rapidly intensified by 10 km downstream. Significant variance also appears at k = 8 and k =16. Melting of deep keels along with ice deformation and artificial smoothing of the ice decreases the channel amplitudes toward the ice front.

The evolution of the cross-shelf thickness variations as the ice moves through the shelf is presented in Figure 8, which shows the Fourier spectrum of cross-shelf variations for each location y along the shelf. Wavenumbers k = 4 and k =12 are dominant at y = 0, as determined by our choice of H0(x), but significant variance develops at wavenumbers k = 8 and k = 16 as melting reshapes the ice. The broadening of the initially sharp spectral peaks at k = 4 and k =12 over the range 0 km < y < 5 km is consistent with the evolution of the basal channels from purely sinusoidal to a more complex wave shape.

Dominant terms in the mass-balance equation

If we decompose the ice velocity, v, as (v − v0) + v0 then, neglecting the surface mass-balance term and the artificial smoothing term (which are small relative to basal melting), the steady-state ice thickness equation, Eqn (6), becomes


We found that in the steady state submarine melting is not balanced primarily by ice deformation (terms involving velocity gradients) but rather by the plug-like advection of ice by the spatially uniform inflow velocity, v0 = 1 km a-1. In other words, the first two terms on the left-hand side of Eqn (21) are found to be rather small, except along the southern edge where velocity gradients tend to thin the ice. In general, the dominant balance is between and − m (with correlation r = 0.99). Since ice deformation is a small effect in the thickness balance, basal channels are essentially the record of spatially variable melting ‘burned’ into the base of the ice as the ice passes like a solid slab over the regions of high melt rates. In this sense, ocean melting acts like a steel- rule die or stencil, imposing a profile on the ice base which is extruded into channels by the rapid seaward motion of the ice mass. A quantitative analysis of this process is given in Appendix B. By y = 10 km, channels are fully formed and subsequently undergo a slow decay due to ice deformation, the gradual melting of the deep keels protruding into warmer ambient water, and artificial smoothing to a small degree. One consequence of this finding is that the existence of basal channels does not imply that the melting itself is distributed spatially in elongated regions. With fast ice flow, it is enough that melting varies significantly in the direction parallel to the grounding line.

In ice shelves with more complex large-scale deforma- tional behavior (e.g. radial spreading), the morphology of channels may be significantly influenced by this large-scale deformation. Our results for Petermann suggest, however, that it is unlikely that velocity gradients at the width scale of the channels themselves provide a significant feedback on channel development or decay.

Origin of high-melt-rate regions

Figure 9 gives an overview of the most important ocean variables for the control case at steady state. Melt is most rapid where the ocean mixed-layer speed is greatest, since for fixed thermal forcing turbulent heat transfer varies linearly with speed. The mixed layer is warmest and saltiest just downstream of the grounding line and deepest at the crests of channels at y ≈ 5 km, before extremely rapid detrainment causes the mixed layer to lose nearly all its volume. Observations (Johnson and others, 2011) confirm that water at mid-depths in the fjord is significantly influenced by under-ice melting, so rapid detrainment rates are plausible. Melt rates correlate closely with ocean mixed-layer speed (r = 0.97) and somewhat less closely (r = 0.90) with thermal forcing, T − T b, since thermal forcing tends to be approximately uniform in the warm deep-water areas that have melt rates m > 20 m a-1.

Fig. 9. Ocean variables for the control case at steady state, (a) Mixed-layer speed, . (b) Mixed-layer thickness, D. (c) Mixed-layer salinity, S. (d) Ocean-induced melt rates, m. (e) Entrainment and detrainment are shown. Positive values show e and negative values show -d. The black contours show the boundary between entraining and detraining regions. (f) Mixed-layer temperature, T. Mixed-layer speed has nearly the same spatial variability as basal melt rate. The ocean mixed layer thickens throughout the eight obvious regions of high melting (south of y = 10 km) before drastically thinning, due to high detrainment rates around y = 8 km. Speed also drops drastically downstream of the high detraining areas as momentum is lost to the ambient ocean. Mixed-layer thickness increases dramatically in the low-melting downstream region (y > 20 km) and tends to flatten the lower surface of the mixed layer (see Fig. 6). Warmer temperatures are associated with higher melting, but temperature alone does not account for the high spatial variability of melt. The most saline mixed-layer water appears at intermediate depths, between y = 3 and y = 15 km, where entrainment of salty water overcomes the freshening effect of melting.

To understand why and where the ocean layer is fast- flowing, we investigate the relative size of terms in the momentum equations, Eqns (9) and (10). In Figure 10 the steady-state values of these terms are shown. In the x – direction there is clearly a dominant geostrophic balance between the Coriolis acceleration and interface flattening terms, giving the approximate relation


Fig. 10. The top row shows terms in the x-transport equation, Eqn (9). (a) Inertial terms, (b) Coriolis acceleration, −DfV. (c) Interfacial-flattening term, g’D(∂a/∂x). d) Detrainment loss, −ḋU. The bottom row shows terms in the y-transport equation, Eqn (10). (e) Inertial terms, (f) Coriolis acceleration, DfU. (g) Interfacial-flattening term, g’D(∂a/∂y). (h) Detrainment loss, −ḋV. The dominant balance in the x-equation is geostrophic (i.e. between the Coriolis acceleration and the interfacial-flattening term, which is the depth-integrated pressure gradient). Inertial terms, eddy diffusion of momentum (not shown), acceleration due to density gradients (not shown), loss of momentum by detrainment and drag forces (not shown) are all small. In the y-equation the balance is between inertial terms, interfacial flattening and detrainment loss. The Coriolis acceleration, eddy diffusion (not shown), density gradient acceleration (not shown) and drag (not shown) are of secondary importance. Axes are distances in km, restricted to the southernmost 20 km to emphasize the region of high melting. Color axes have units m2 s-2.

In the y –direction the balance is different, with three dominant terms yielding


According to Eqn (23), steep along-channel slopes in the ambient interface, z = a(x, y), accelerate water northward up the channels to y « 6 km, where intense melting causes the mixed layer to detrain rapidly, at which point the mixed- layer volume, and hence momentum, is lost into the ambient ocean. Lacking a 3-D ocean we are subsequently unable to account for this momentum. Relation (22) states that as large velocities, V, develop, a tilt in the ambient interface develops to geostrophically balance the along-channel jets. Indeed, in Figure 6 we see that in the channel crests at y = 3, 5 and 7 km the ambient interface (red curves) has a pronounced tilt just below the crests. It is this banking of the mixed- layer water against the eastern side of the channel crests that supports the along-channel jets responsible for rapid melting (Fig. 5a and b). The melt rates of ~30 m a-1 occurring outside of these jets and south of y = 10 km are associated with eastward flow cutting diagonally across isobaths. The mixed- layer thickness, D, is rather small in these regions and, since terms in the transport equation are weighted by D, Eqns (22) and (23) do not explain the dynamics in these regions.

Partial isolation of circulation in a channel

The circulation within a given channel is not independent of the neighbouring channels. In Figure 11 cross-shelf water transport, DU, is plotted against the ice base gradient, ∂b/∂x, for all locations on the flanks or spine of an ice keel (locations where 2b/∂x2 > 0). The clear positive correlation shows that water tends to be transported in the direction of the ice base’s gradient. Note, however, that points where ∂b/∂x = 0 tend to have DU > 0. Figure 11 therefore shows that significant amounts of water are transported eastward across the deep ice keels, so that circulation in an individual channel is not isolated from flow in other channels.

Fig. 11. Cross-shelf transport plotted against cross-shelf gradient, for locations on the flanks and base of keels (∂2b/∂x2 > 0).

Perturbation experiments

Many experiments involving perturbations to the control case steady state were performed. The sets of experiments we present in detail involved perturbations to the thickness profile, H0(x), of ice entering along the southern edge, perturbations to the temperature, Ta(z), of the ambient water and perturbations to the amount of subglacial discharge, s, entering the ocean near the southern boundary.

In general, scenarios where the ocean melts completely through the shelf are beyond the capability of our model, and indeed cases where melting creates exceptionally deep and steep channels are also impossible to simulate, since in such cases we are unable to find converged ice velocities. We present only experiments where converged results were obtained and no other model failures occurred.

Thickness perturbations

First, we show the results of modifying the inflowing ice thickness profile on the southern edge from the ‘4 plus 12’ profile in the control case to a case where the inflow profile is sinusoidal with single wavenumber k and amplitude A. We investigated values k ∊ 0,1,2, …,15 and A ∊ 25.0,50.0 m. The results of these experiments appear in Figure 12, in which the area-integrated melting, expressed as a percentage of the total volumetric flux across y = 0km (12 km3 a-1 in all cases), is plotted against the wavenumber, k, of the perturbation. Note that surface ablation of 1 ma-1 accounts for 7% of the total mass loss already, so total melting over ~90% makes the ice shelf unviable.

Fig. 12. The blue solid line (25 m perturbations) and blue crosses (50m perturbations) give the area-integrated applied melt rates as a percentage of the total volumetric flux of ice across the inflow edge (12 km3 a-1 for all cases). The horizontal axis shows the wavenumber of the sinusoidal perturbation imposed on the ice along the southern boundary. The black solid line (25 m perturbations) and black crosses (50m perturbations) show the time derivative of total ice-shelf volume (again as a percentage of 12 km3 a-1) when the simulation stopped, either because a steady state was reached (all k > 5) or because the model could not proceed (all k < 5). For k < 5 the homotopy time interval, tf, had not yet passed (see Table 1) and for these cases the red solid line (25 m) and red crosses (50 m) show the area integral of m, while the blue line and crosses show the area integral of 7/11 at the end of the simulation.

We draw two conclusions from these experiments. First, the amplitude of the initial perturbation has very little effect on the eventual total melt rate, at least in cases where a steady state is achieved. Secondly, higher-wavenumber perturbations lead to decreased overall melting, and wavenumbers k < 4 generally lead to so much mass loss (exceeding 100% of the influx rate, if surface ablation is included) that the ice becomes excessively thin either along the eastern boundary or else along channel crests, leading to a breakdown of the model. The time of model failure for these cases is given in Table 1, along with the spin-up time required (including the 100 year homotopy interval; see Appendix A) for other cases to reach a steady state.

To the extent to which our simplified model represents the real ice/ocean system, this suggests that Petermann’s ice shelf might experience much more melting if irregularities in thickness were absent at the grounding line. If Petermann is currently stable, as we have assumed, the deep channels we observe may be a necessary defense against excessive melting, making the present ice shelf’s areal extent viable. We note also that a wavenumber of 6 or 7 (wavelength of ~3 km) yields an overall melt rate of 80% of mass influx, which is approximately the correct amount of melting for Petermann Gletscher (Johnson and others, 2011).

Temperature perturbations

Recent warming in nearby Nares Strait (Munchow and others, 2011) and the likelihood that warming of coastal waters has contributed to ice retreat and acceleration elsewhere around Greenland (Holland and others, 2008; Murray and others, 2010; Christoffersen and others, 2011; Straneo and others, 2011) motivated a set of experiments in which the ambient water properties, Ta(z) and Sa(z) were changed.

In one set of experiments the temperature of the surface water was changed by adding an increment, AT, to the uppermost two control points in the temperature profile shown in Figure 4, affecting all ambient temperatures above 1 75 m depth. A warming of 0.3°C caused already thin ice near the ice front to further thin by tens of meters (Fig. 13), but thicker ice in contact with deeper water was generally unaffected. Because thinning occurred mostly near the northern outflow edge, the outgoing ice flux decreased significantly, such that melting accounted for an additional 6.6% of total mass loss. Note that the total melting, total surface ablation (which is fixed) and the flux of ice across the ice front must all sum to 12 km3 a-1, so an increase in total melting is always indicated by thinning of the ice exiting along the northern edge of the shelf. Thinning of ice at intermediate depths alone indicates a rearrangement in the pattern of melting, but not necessarily an increase in total melting.

Fig. 13. Four temperature perturbation experiments were performed. In each panel the perturbed ice drafts are plotted against the initial ice drafts. Points falling on the black lines indicate locations where the ice draft was unchanged. (a) The deep ambient water was warmed by 0. 3ºC. We see significant thinning of the ice at depths from 450 to 100m. (b) Warming the upper ambient water column by the same amount leads to significant thinning of the shallowest ice and hence to a large increase in the total amount of melting (an extra 6.6% of total mass loss). Deeper ice is little affected by the perturbation. (c) Cooling the deep water by 0. 1 C leads to ice thickening and an overall melt decrease, which is linearly consistent with the 0. 3º C warming response shown in (a). (d) Cooling by 0. 3ºC, however, leads to a surprising result where, although some ice thickens at intermediate depths, ice in other areas thins dramatically, leading to model failure before a steady state is reached. At this time the total melting exceeded the control case melting by 0.5 km3 a 1.

Next, perturbations to the subsurface water (below 200 m depth) were carried out. The greatest subsurface warming we were able to test without encountering a model failure from excessive melting was AT = 0.3oC, which led to significant thinning of the ice at intermediate depths but relatively little thinning of the thinnest ice (Fig. 13a). The outgoing flux of ice across the northern edge therefore decreased by relatively little and so the overall amount of melting increased by a correspondingly small amount, from 78.9% of the total rate of mass loss to 80.4%. The main effect of subsurface warming of this degree is that the area of rapid melting shifts southward, meaning the ice base becomes steeper and the base retreats into cold surface water at a more southerly location.

Cooling the deep water by 0.1oC led to thickening of the ice at intermediate depths and a decrease in overall melting from 78.9% of total mass loss to 78.3%, which is linearly consistent with the response to a 0. 3oC subsurface warming. Cooling the subsurface water by 0.3oC led, surprisingly, to a reorganization of the melt pattern, such that melting in some areas became large enough to melt nearly through the ice 12 km from the grounding line. An instantaneous view (Fig. 14) of the ice draft and melt rates at the time of model failure shows that four regions of extremely high melting formed between y = 5 and y = 10 km. In the control case, deepening of the k = 12 undulations serves to reduce the speed of the ocean layer and reduces the amount of melting that would occur in the otherwise- unviable k = 4 undulations. It appears that the colder water in this experiment failed to deepen the k = 12 undulations enough to prevent the development of high melting at the k = 4 scale. In this way, the cooling of deep water pushed the system towards the unviable k = 4 case (Fig. 12). This is a very specific result which is unlikely to directly apply to Petermann Gletscher, but it shows that the measure of defense provided by narrow channels can be overwhelmed by an initially slight reorganization of the melt-rate pattern due to an otherwise innocuous temperature change.

Fig. 14. (a) Colors show ice draft at time of model failure during the adjustment to a 0.3°C subsurface cooling in the ambient layer. Overlaid arrows are the mixed-layer velocity. (b) Melt rates at time of model failure. (c) Channel depths, again at the time of model failure. Note the four regions beginning around y = 12 km, where high melting has created deeply incised channels. (a) shows that the draft in these regions is < 50 m.

Subglacial discharge experiments

Several experiments were performed in which positive subglacial discharge rates, ṡ, were used. Beginning with the control case steady state, where = 0 everywhere, we increased ṡ uniformly in a rectangular region south of y = 1. 0 km, such that the area integral of ṡ took the value 0.5, 1.0, 1.5 or 2.0 km3 a 1. In general, increasing total discharge from 0 to 1.5 km3 a 1 increases the total amount of melting in the resulting steady state (see Fig. 15). When total discharge exceeds 1. 5km3 a 1, instantaneous melt rates are even larger and excessive thinning of the ice prevents the model from reaching a steady state. As with subsurface warming, increasing subglacial discharge does not cause the shallowest ice (above 50m depth) to thin very much. There is, however, significant thinning at intermediate depths, where the ice base rises by a few tens of meters for a discharge rate of 0.5 km3 a 1 and up to 100m of thinning for a discharge rate of 1.5 km3 a 1. In these cases the channeldeepening region simply shifts southward with increasing rates of discharge, which explains why intermediate depths are most affected.

Fig. 15. Increase in total steady-state melting in response to the addition of subglacial discharge. The horizontal axis is the area- integrated d and the vertical axis is the additional area-integrated m relative to the control case, for which ṡ = 0.

Multiple steady states

Given a set of boundary conditions and forcing terms it is possible that the steady state achieved by the model depends on the initial state of the model. We investigated this possibility in a limited way and found that, at least for subglacial discharge perturbations, the final steady state is sensitive to the initial ice state. In particular, if subglacial discharge is activated beginning with ice in the pre-control state, the eventual steady-state melt rates and ice thickness (Fig. 16, first column) differs from that of a simulation in which discharge is activated with ice in the control case steady state (Fig. 16, third column). Figure 16 demonstrates that more than one steady state exists for this set of boundary conditions and external forcing. We probed the stability of the pre-control steady state by applying the same amount of subglacial discharge reorganized into four jets by weighting s by the function . We then recover a steady state similar to the control case regime (Fig. 16, second column). This experiment shows that the specific distribution of subglacial discharge is enough to nudge the system from one stable regime to another.

Fig. 16. (a) Steady-state melt patterns and (d) ice thickness when subglacial discharge totaling 1.5 km3 a-1 is added uniformly to the region y < 1 km, beginning with the pre-control state. If the same discharge is applied beginning with the control state the resulting steady-state melt pattern (c) and ice thickness (f) differ from (a) and (d), although the forcing and boundary conditions are the same. If the same volume of discharge is organized into four broad jets and applied to the pre-control case, the eventual steady-state melt rate (b) and thickness patterns (e) exhibit the dominant k = 4 character, as exhibited by (c, f).

It appears that the origin of multiple steady states is the nonlinearity of the ocean circulation which creates a competition between amplification of basal undulations on different width scales. When subglacial discharge is applied to the pre-control case, it induces enough melting at the k =12 scale that the resulting channels are able to control the flow of water and prevent amplification of the k = 4 channels. However, when subglacial discharge is applied to the control case, deep undulations at wavenumber k = 4 already exist and the enhanced melting is insufficient to allow the k =12 mode to control the flow of water. When subglacial discharge is organized into four jets, the k = 4 mode grows quickly enough to suppress the enhancement of the k =12 mode, even though the system begins in the pre-control state.

Our use of the homotopy method (Appendix A), which tends to break the physical connection between the initial state and final state, certainly casts doubt on this explanation, yet the existence of distinct steady states is not at all invalidated by the use of the homotopy method, since the results shown all correspond to times well after Y =1.0 was attained.

Entrainment perturbations

The entrainment scheme used by the ocean model is extremely important in determining the heat flux to the ice and, consequently, the melt rates and ice geometry. Perturbations of the parameter, 1, appearing in the TKE source term in Eqn (11) led to very large adjustments in the ice geometry. Changing ¡1 from 2.5 (the value used in the control case) to 1.0 deepened and steepened basal channels to the point of causing model failure. However, adjusting ¡1 to 5.0 led to a steady state where even the incipient channels in H0(x) were erased by y = 15 km. This occurred because enhanced entrainment deepened the ocean mixed layer enough to override the steering influence of the incipient channels. Broad swaths of melting then appeared that eventually removed any ice base irregularities further downstream. It may be that the real system is correspondingly sensitive to changes in ocean turbulence, or it may be that our model only gives results that appear correct in a narrow parameter space in which errors cancel in a fortuitous manner.

All the results described here used a Niiler-Kraus-type scheme for TKE dissipation, but we also experimented with Gaspar’s ‘CMO’ scheme (both schemes are described by Gaspar, 1988). We were able to obtain qualitatively similar results using the CMO scheme, which we do not discuss further here. A third scheme, from Zilitinkevich and Mironov (1996), was also used, which, instead of diagnosing an entrainment rate, diagnoses an ‘equilibrium mixed-layer depth’ directly. In this approach the entrainment or detrainment rate is simply determined by the requirement to relax the mixed-layer depth to the new diagnosed depth in the current time-step. Again, using fairly liberal tuning of the non-dimensional parameters in their scheme we were able to obtain a coupled steady state qualitatively similar to the control case presented above.

In contrast, the entrainment scheme attributed to Kochergin by Holland and Feltham (2006), which does not permit detrainment, did not produce channelized topography in the parameter space we explored. It generally produced too much melting along the eastern boundary, where Coriolis acceleration tends to pile up mixed-layer water. For schemes permitting detrainment, melting along the eastern wall causes the mixed layer to detrain, moving the accumulated warm water away from contact with the ice base and therefore preventing excessive melting. Since observations also suggest that melt-influenced water appears at intermediate depths outside the fjord mouth (Johnson and others, 2011), we conclude that detrainment is probably significant in the circulation under the Petermann Ice Shelf.

Other perturbation experiments

In addition to the experiments described above in detail, we briefly report the results of other experiments. Unless stated otherwise, we report results where the adjustment to a perturbation reached a final steady state.

Increasing turbulent horizontal mixing by increasing Ah and Kh to 100 m2s-1 shuts off channel formation by smoothing and broadening the regions of intense melting. Increasing the mixed-layer minimum thickness, D0, from 1.0m to 10.0m also suppresses channel formation by broadening the regions of rapid melt so significant melting occurs on the keels.

When U*0, the contribution to the friction velocity from unmodeled flow, is decreased from 2.5 x 10-3 cms-1 to 1.0 x 10-3 cms-1, the total amount of melting decreases to 71% of the influx rate (from 78.9% in the control case) but the maximum channel depth remains at ~250m. However, doubling U*0 to 5.0 x 10-3cms-1 causes channels to be eroded completely and the ice is melted through along the eastern boundary before steady state can be achieved.

Speeding up the influx of ice across the southern edge causes the regions of intense melting to weaken and elongate in the y-direction. If v0 is increased to 1500 m a-1, channels become shallower, with a maximum keel-to-crest amplitude of ~175m, and, although the total amount of melting increases, the amount of melting relative to the influx of ice decreases.

Some parameter perturbations had very little influence on the amount and pattern of melting. Increasing the lateral stress T 0 from 25 kPa to 50kPa decreases the maximum ice speed from ~1100ma-1 to ~1050ma-1, and differences in ice geometry and melt rates are extremely minor, with channel deepening occurring perhaps 1 km further downstream. Similarly, warming the ice interior by 5oC has the effect of steepening the ice very slightly just north of the sourthern boundary, which compresses the melt regions and shifts the fully channelized part of the ice by ~1 km up- glacier. Cooling by the same amount has a slightly smaller effect with the opposite sign.

We found that quite specific ranges of values were necessary for certain key parameters (i.e. Ah Kh U*0 and j) in order for the model to qualitatively match the observed geometry of Petermann Gletscher. Only Ah, and Kh, appear in other applications of the same ocean model. Payne and others (2007) and Holland and others (2009) used Ah = Kh, = 100 m2 s-1 with a grid spacing of 1 km. Holland and others (2007) used Ah, = Kh, = 500 m2s-1 with a grid spacing of 2 km. These are all considerably larger than our value of Ah , = Kh, = 10 m2 s-1 (with grid spacing of 250 m), even if scaled by the grid spacing. Again, the narrowness of the parameter space in which our model exhibits realistic behavior may indicate fortuitous model error cancellations, or it may reflect a true sensitivity in the real system. The uniqueness of Petermann within Greenland perhaps supports the latter. Only improved observations of sub-shelf ocean turbulence can improve confidence in models like ours.


To summarize our results, we have found that ocean- induced melting, combined with rapid ice flow, is enough to produce the channelized basal morphology observed at Petermann Gletscher and other glaciers. The presence of narrow channels likely reduces the total amount of melting by suppressing melting at broader spatial scales. Ice deformation is not involved in the dynamical balance responsible for channel deepening.

We were not able to determine any dynamical process that ultimately selects the width scale of channels, although it appears that irregularities imprinted on the ice by the bedrock are preferentially amplified by ocean melting. Calculation of , the Rossby radius for the ocean mixed-layer model, gives a range of values between 0.1 and 2.7km for the control case steady state. In the regions of high melt the Rossby radius is ~2.5 km, which is consistent with the width scale of the fast-flowing jets responsible for the high melt. Although there is agreement in scales here, no mechanism relating the Rossby length to the width of channels is immediately apparent.

The ice geometry is sensitive to ocean temperatures, although catastrophic thinning of the ice seems more likely in response to surface warming of the ocean than to subsurface warming. Surface warming of a few tenths of a degree in our model leads to melting that penetrates the ice tens of kilometers inland of the ice front. This result is questionable, however, given our incorrect salinity profile for the ambient ocean and our inability to represent the intrusion of cold surface water under the ice shelf. Modeling the fully 3-D ocean is required to address this. Warming of the subsurface waters by a few tenths of a degree does not lead to melting that penetrates the ice shelf anywhere near the grounding line, though some thinning of the already thin seaward ice does occur.

It has been suggested (Rignot and Steffen, 2008) that basal channels may make Petermann and other similarly channelized ice shelves more susceptible to disintegration. The present model cannot adequately test this hypothesis, since we treat ice as a viscous material and neglect the brittle character of ice which is necessary for sudden mechanical failure. To project, for instance, the fate of the Pine Island Ice Shelf as the grounded glacier continues its well-observed (Rignot, 2008; Pritchard and others, 2009; Joughin and others, 2010) thinning and acceleration, numerical models will need to account for the brittle fracture-permitting rheology of ice, as well as smooth ice morphology resulting from ice/ocean interaction. In fact, the two problems are related, since patterns of surface and basal crevasses have been related to the locations of basal channels under the Pine Island Ice Shelf (Vaughan and others, 2012). Stresses tending to induce such crevasses are non-hydrostatic and therefore cannot be studied with the present HOM model.

The discharge of fresh water along the grounding line is found to have a significant effect on submarine melting in the model, with a total discharge of 2.0 km3 a-1, deepening basal channels enough to cause model failure. Lesser amounts of subglacial discharge distributed uniformly led to steady states with enhanced channelization. Multiple steady states are possible in the coupled model and we found that concentrating subglacial discharge into four broad jets is enough to nudge the system from one steady state (which lacks channels of wavenumber k = 4) to another featuring strong variance at wavenumber k = 4.

Excessive turbulent diffusion in the ocean (large values of A, and K,) is enough to shut off channel formation, as is the variation of certain entrainment parameters (j and U*0) that are otherwise poorly constrained. There is therefore an acute need for better observations of turbulence in sub- ice-shelf environments if modelers are to have a chance of correctly simulating long-term ice-shelf melting in warm ocean conditions. A perhaps deeper concern regarding these entrainment schemes is their universal assumption that the TKE budget is closed within each vertical column. On the small grids we have found necessary to resolve complex ice geometry, ocean currents may traverse a gridcell in 15 min, making the assumption of a vertically balanced TKE budget at least questionable.

For ice modelers, we emphasize the need for models to cope with extreme variability in ice thickness. For example, it may be necessary to permit ice to vanish completely from interior cells (i.e. cells potentially far from the ice front) if models are to be robust enough to simulate extreme retreat of marine glaciers. Finally, with melting varying on small spatial scales it seems likely that the adjustment to local hydrostatic equilibrium occurs on a timescale that is comparable to, or larger than, the residence time of the ice in the ice shelf. The use of the Blatter-Pattyn equations for rapidly melting ice shelves may lead to significant errors in ice geometry and stress state.


We thank Andrew Fleming and Tom Milgate for providing Figure 1. Bill Lipscomb provided valuable guidance regarding model coupling. Adrian Jenkins made helpful suggestions regarding entrainment and subglacial discharge. This study was funded by the IMPACTS project sponsored by the US Department of Energy’s Program in Biological and Environmental Research.


Baines, PG (2008) Mixing in downslope flows in the ocean – plumes versus gravity currents. Atmos.-Ocean, 46(4), 405419 (doi: 10.3137/ao.460402)
Bindschadler, R, Vaughan, DG and Vornberger, P (2011) Variability of basal melt beneath the Pine Island Glacier ice shelf, West Antarctica. J. Glaciol., 57(204), 581595 (doi: 10.3189/002214311797409802)
Blatter, H (1995) Velocity and stress fields in grounded glaciers: a simple algorithm for including deviatoric stress gradients. J. Glaciol., 41(138), 333344
Christoffersen, P and 7 others (2011) Warming of waters in an East Greenland fjord prior to glacier retreat: mechanisms and connection to large-scale atmospheric conditions. Cryosphere, 5(3), 701714 (doi: 10.5194/tc-5–701–2011)
Cuffey, KM and Paterson, WSB (2010) The physics of glaciers, 4th edn. Butterworth-Heinemann, Oxford
Dukowicz, JK, Price, SF and Lipscomb, WH (2010) Consistent approximations and boundary conditions for ice-sheet dynamics from a principle of least action. J. Glaciol., 56(197), 480496 (doi: 10.3189/002214310792447851)
Durran, DR (1999) Numerical methods for wave equations in geophysical fluid dynamics. Springer, New York
Fricker, HA, Coleman, R, Padman, L, Scambos, TA, Bohlander, J and Brunt, KM (2009) Mapping the grounding zone of the Amery Ice Shelf, East Antarctica using InSAR, MODIS and ICESat. Antarct. Sci., 21(5), 515532 (doi: 10.1017/S095410200999023X)
Gaspar, P (1988) Modeling the seasonal cycle of the upper ocean. J. Phys. Oceanogr., 18(2), 161180 (doi: 10.1175/1520- 0485(1988)018<0161:MTSOTT >2.0.CO;2)
Grosfeld, K and Sandhager, H (2004) The evolution of a coupled ice shelf-ocean system under different climate states. Global Planet. Change, 42(1–4), 107132 (doi: 10.1016/j.gloplacha.2003.11.004)
Hansen, MO, Nielsen, TG, Stedmon, CA and Munk, P (2012) Oceanographic regime shift during 1997 in Disko Bay, Western Greenland. Limnol. Oceanogr., 57(2), 634644 (doi: 10.4319/lo.2012.57.2.0634)
Heimbach, P and Losch, M (2012) Adjoint sensitivities of sub-ice- shelf melt rates to ocean circulation under the Pine Island Ice Shelf, West Antarctica. Ann. Glaciol., 53(60 Pt 1), 5969 (doi: 10.3189/2012AoG60A025)
Holland, PR and Feltham, DL (2006) The effects of rotation and ice shelf topography on frazil-laden ice shelf water plumes. J. Phys. Oceanogr., 36(12), 23122327 (doi: 10.1175/JPO2970.1)
Holland, DM and Jenkins, A (2001) Adaptation of an isopycnic coordinate ocean model for the study of circulation beneath ice shelves. Mon. Weather Rev., 129(8), 19051927
Holland, DM, Jacobs, SS and Jenkins, A (2003) Modelling the ocean circulation beneath the Ross Ice Shelf. Antarct. Sci., 15(1), 1323 (doi: 10.1017/S0954102003001019)
Holland, PR, Feltham, DL and Jenkins, A (2007) Ice shelf water plume flow beneath Filchner-Ronne Ice Shelf, Antarctica. J. Geophys. Res., 112(C5), C05044 (doi: 10.1029/2006JC003915)
Holland, DM, Thomas, RH, de Young, B, Ribergaard, MH and Lyberth, B (2008) Acceleration of Jakobshavn Isbræ triggered by warm subsurface ocean waters. Nature Geosci., 1(10), 659664 (doi: 10.1038/ngeo316)
Holland, PR, Corr, HFJ, Vaughan, DG, Jenkins, A and Skvarca, P (2009) Marine ice in Larsen Ice Shelf. Geophys. Res. Lett., 36(11), L11604 (doi: 10.1029/2009GL038162)
Jenkins, A and Holland, DM (2002) A model study of ocean circulation beneath Filchner-Ronne Ice Shelf, Antarctica: implications for bottom water formation. Geophys. Res. Lett., 29(8), 1193 (doi: 10.1029/2001GL014589)
Jenkins, A and 6 others (2010) Observations beneath Pine Island Glacier in West Antarctica and implications for its retreat. Nature Geosci., 3(7), 468472 (doi: 10.1038/ngeo890)
Johnson, HL, Münchow, A, Falkner, KK and Melling, H (2011) Ocean circulation and properties in Petermann Fjord, Greenland. J. Geophys. Res., 116(C1), C01003 (doi: 10.1029/2010JC006519)
Joughin, I and Alley, RB (2011) Stability of the West Antarctic ice sheet in a warming world. Nature Geosci., 4(8), 506513 (doi: 10.1038/ngeo1194)
Joughin, I, Smith, BE and Holland, DM (2010) Sensitivity of 21st century sea level to ocean-induced thinning of Pine Island Glacier, Antarctica. Geophys. Res. Lett., 37(20), L20502 (doi: 10.1029/2010GL044819)
Lemieux, J-F and 6 others (2011) Implementation of the Jacobian- free Newton-Krylov method for solving the first-order ice sheet momentum balance. J. Comput. Phys., 230(17), 65316545 (doi: 10.1016/
Lipscomb, WH and Hunke, EC (2004) Modeling sea ice transport using incremental remapping. Mon. Weather Rev., 132(6), 13411354 (doi: 10.1175/1520–0493(2004)132<1341: MSITUI>2.0.C0;2)
MacAyeal, DR, Rommelaere, V, Huybrechts, P, Hulbe, CL, Determann, J and Ritz, C (1996) An ice-shelf model test based on the Ross Ice Shelf, Antarctica. Ann. Glaciol., 23, 4651
Mankoff, KD, Jacobs, SS, Tulaczyk, SM and Stammerjohn, SE (2012) The role of Pine Island Glacier ice shelf basal channels in deep-water upwelling, polynyas and ocean circulation in Pine Island Bay, Antarctica. Ann. Glaciol., 53(60 Pt 1), 123128 (doi: 10.3189/2012AoG60A062)
Motyka, RJ, Truffer, M, Fahnestock, M, Mortensen, J, Rysgaard, S and Howat, I (2011) Submarine melting of the 1985 Jakobshavn Isbr* floating tongue and the triggering of the current retreat. J. Geophys. Res., 116(F1), F01007 (doi: 10.1029/2009JF001632)
Mueller, RD, Padman, L, Dinniman, MS, Erofeeva, SY, Fricker, HA and King, MA (2012) Impact of tide-topography interactions on basal melting of Larsen C Ice Shelf, Antarctica. J. Geophys. Res., 117(C5), C05005 (doi: 10.1029/2011JC007263)
Munchow, A, Falkner, KK, Melling, H, Rabe, B and Johnson, HL (2011) Ocean warming of Nares Strait bottom waters off northwest Greenland, 2003–2009. Oceanography, 24(3), 114123 (doi: 10.5670/oceanog.2011.62)
Murray, T and 10 others (2010) Ocean regulation hypothesis for glacier dynamics in southeast Greenland and implications for ice sheet mass changes. J. Geophys. Res., 115(F3), F03026 (doi: 10.1029/2009JF001522)
Pattyn, F (2003) A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res., 108(B8), 2382 (doi: 10.1029/2002JB002329)
Payne, AJ and Dongelmans, PW (1997) Self-organization in the thermomechanical flow of ice sheets. J. Geophys. Res., 102(B6), 12 21912 233 (doi: 10.1029/97JB00513)
Payne, AJ, Holland, PR, Shepherd, AP, Rutt, IC, Jenkins, A and Joughin, I (2007) Numerical modeling of ocean-ice interactions under Pine Island Bay’s ice shelf. J. Geophys. Res., 112(C10), C10019 (doi: 10.1029/2006JC003733)
Price, SF, Payne, AJ, Howat, IM and Smith, BE (2011) Committed sea-level rise for the next century from Greenland ice sheet dynamics during the past decade. Proc. Natl Acad. Sci. USA (PNAS), 108(22), 89788983 (doi: 10.1073/pnas.1017313108)
Pritchard, HD, Arthern, RJ, Vaughan, DG and Edwards, LA (2009) Extensive dynamic thinning on the margins of the Greenland and Antarctic ice sheets. Nature, 461(7266), 971975 (doi: 10.1038/nature08471)
Pritchard, HD, Ligtenberg, SRM, Fricker, HA, Vaughan, DG, Van den Broeke, MR and Padman, L (2012) Antarctic ice-sheet loss driven by basal melting of ice shelves. Nature, 484(7395), 502505 (doi: 10.1038/nature10968)
Rignot, E (2008) Changes in West Antarctic ice stream dynamics observed with ALOS PALSAR data. Geophys. Res. Lett, 35(12), L12505 (doi: 10.1029/2008GL033365)
Rignot, E and Kanagaratnam, P (2006) Changes in the velocity structure of the Greenland Ice Sheet. Science, 311(5673), 986990 (doi: 10.1126/science.1121381)
Rignot, E and Steffen, K (2008) Channelized bottom melting and stability of floating ice shelves. Geophys. Res. Lett, 35(2), L02503 (doi: 10.1029/2007GL031765)
Rutt, IC, Hagdorn, M, Hulton, NRJ and Payne, AJ (2009) The Glimmer community ice sheet model. J. Geophys. Res, 114(F2), F02004 (doi: 10.1029/2008JF001015)
Schoof, C and Hindmarsh, RCA (2010) Thin-film flows with wall slip: an asymptotic analysis of higher order glacier flow models. Q. J. Mech. Appl. Math, 63(1), 73114 (doi: 10.1093/qjmam/hbp025)
Stewart, C, Rignot, E, Steffen, K, Cullen, K and Ruff, R (2004) Basal topography and thinning rates of Petermann Gletscher, northern Greenland, measured by ground-based phase-sensitive radar. FRISP Rep. 15
Straneo, F and 6 others (2011) Impact of fjord dynamics and glacial runoff on the circulation near Helheim Glacier. Nature Geosci., 4(5), 322327 (doi: 10.1038/ngeo110)
Thoma, M, Jenkins, A, Holland, D and Jacobs, S (2008) Modelling circumpolar deep water intrusions on the Amundsen Sea continental shelf, Antarctica. Geophys. Res. Lett, 35(18), L18602 (doi: 10.1029/2008GL034939)
Thomas, R, Frederick, E, Krabill, W, Manizade, S and Martin, C (2009) Recent changes on Greenland outlet glaciers. J. Glaciol., 55(189), 147162 (doi: 10.3189/002214309788608958)
Vaughan, DG and 8 others (2012) Subglacial melt channels and fracture in the floating part of Pine Island Glacier, Antarctica. J. Geophys. Res., 117(F3), F03012 (doi: 10.1029/2012JF002360)
Zilitinkevich, SS and Mironov, D (1996) A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer. Bound.-Layer Meteorol., 81(3–4), 325351 (doi: 10.1007/BF02430334)

Appendix A: Numerical Robustness Considerations

Since we are investigating ice-shelf failure from excessive melting, we regularly encounter model failures and yet we would like to be able to interpret model failure as physically meaningful in some cases. However, in certain model configurations failures occur which are essentially artificial. For instance, if we spin-up the model beginning with ice having the profile of an unconfined shelf and floating on top of warm water, then a region of thin ice develops just seaward of the grounding line, due to high melt rates there. Because the initially thick ice has not yet entirely advected out of the domain, a deep indentation in the ice base results, which can lead to divergence in the ice velocity solver routine. Or, in some perturbation experiments, transient melting occurs that is rapid enough to reduce the ice to just a few meters thickness either along the eastern boundary or on the eastern sides of channel crests. This induces model failure, yet does not allow us to rule out the existence of a steady state for the given boundary conditions and forcing. In such cases we found that a simple homotopy method improved the robustness of numerical convergence in the ice velocity solver throughout the transient adjustment. The method is simply to apply a melt rate, 71V, to the ice instead of the rate, m, determined by the ocean model, where 0 < 7 < 1 is gradually made to approach 1 over a fixed time interval. In this way the ice is allowed to gradually adapt to large melt rates without developing artificially large thickness gradients which induce divergence in the ice velocity solver.

In general, we let γ depend on the elapsed time, t, and the local ice draft, d, according to

where H (x) is the Heaviside step function and

is a minimum melt depth threshold. By defining γ in this way, γ(t, d) = 0 whenever the ice draft, d, is less than the cut-off value, dmin(t), which decreases from dmin(0) = d0 to 0.0 m over a time tf. In this way, melting is shut off for ice that is too thin at a given time, but thinner and thinner ice experiences some melting as time proceeds. Also, for d > dmin(t), γ varies linearly from γ(0, d) = f0 to γ (tf, d) = 1.0 over the same time period, tf. When t reaches tf the ice feels the full melt rate, m, so any steady state achieved after tis time is a true solution. Forthe main spin-up we used d0 = 0, f0 = 0.75 and tf = 100 years. For the thickness perturbation experiments we used d0 = 50.0 m, f0 = 0.75 and tf = 100 years. Using this homotopy method, the intermediate states in the transient solution are therefore partially unphysical, yet the resulting steady state (after t > tf) is a true steady solution of the coupled equations.

The diffusion term, k(∂H/∂x) appearing in Eqn (6) is, similarly, a numerical device that we found was necessary to prevent the ice from developing steep gradients and cusps that sporadically prevent the velocity solver from converging. It is not meant to capture any known rheological properties of ice, although the behavior of ice when it develops extremely steep thickness variations needs to be further investigated. A value of k = 5.0 x 103m2a-1 was found to be large enough to suppress sharp gradients at the gridscale without suppressing the formation of channels at the resolved scale. If k is smaller than x 103m2a-1 we find intense melting usually leads to steep ice gradients that prevent the convergence of the Picard iteration. A term of the form k(∂2H/∂y2) could also have been added to make the artificial smoothing isotropic but this would contribute only a few centimeters per year of ice thickness change. In general we found that ice thickness is influenced much more by melting than by this artificial smoothing. Furthermore, locations with the largest melt rates are generally being thickened by artificial smoothing which means that channelization would be somewhat more pronounced in the absence of smoothing. However, we find that the term Eqn (6) is typically of the same size as the artificial smoothing term (a few meters per year) and both terms tend to fill in channels, which means that we are not able to isolate the effect of transverse ice deformation on channel development.

Appendix B: Channel Formation

Here we present a simple analysis explaining how the balance v0(∂H/∂y)—ṁ gives rise to channels in our control case. Channels are simply the presence of significant basal slopes in the x-direction which persist for a considerable distance in the y-direction. Stated quantitatively, we distinguish channels if ∂b/∂x takes on large values while L(∂2b/∂y∂x), the rate at which those slopes change along the flow direction, remains small, for some length scale, L, corresponding to the length of the channels. Here b = − (ρIo)H is the elevation of the ice base, as in Figure 2. First, we may write formally

where x1 is at the crest of a channel and x0 is at a neighboring keel, while y1 is a position down-glacier relative to y0. For convenience, denote the channel amplitude, b (x1, y) − b (x0, y), by C (y). Then if v0(∂H/∂y) ≈−ṁ we have


In other words, channels deepen by the cumulative effect of melt rates which differ in the cross-shelf direction.

In our model simulations, we typically have ṁ(x1, y) −ṁ (x0, y) ≈ 50 m a-1 near the southern inflow edge and this difference persists for ~5 km along the direction of flow, so we ought to have 5 km = 220m. This is, indeed, about the amount of channel deepening the model simulates and, given the 50 m inflow undulations, this yields a maximal channel wall gradient of |(∂b/∂x)| ≈ 0.18 if we take the keel-to-crest distance Δx = x1 − x0 to be 1500m. The steepest channel wall gradients at the Petermann Ice Shelf are |(∂b/∂x)| = tan θ ≈ 0.36, for a channel side-wall inclination angle of d = 20°. This corresponds to a channel ‘opening angle’ of 140°, as reported by Rignot and Steffen (2008), which is the angle with vertex at the crest of a channel and legs aligned with the channel side-walls. Therefore the steepest channel sidewalls at the Petermann Ice Shelf are still about twice as steep as the model’s channel side-walls.

Equation (B1) then gives

By y = 10 km, melt rates differ by at most 15 m a-1 over a distance of Δx = 2 km so if L = 5 km we have . Therefore, the advective balance, implies that, beyond y = 10 km, the channel shape should not change very much over any 5 km stretch, in the direction of flow, which is indeed accurate for the control case (Fig. 5).