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 (Reference Payne and DongelmansPayne and Dongelmans, 1997) to a uniform ice temperature, T₁ , which is held constant. In the control run we set T₁ = −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 Reference Dukowicz, Price and LipscombDukowicz and others (2010) and Reference Schoof and HindmarshSchoof 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 m² 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 p₁ = 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 Reference Lipscomb and HunkeLipscomb 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 = −H (ρ₁/ρ₀), with ρ₀ = 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 (Reference Johnson, Münchow, Falkner and MellingJohnson 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/η)τ₀, where τ₀ 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 (Reference Cuffey and PatersonCuffey 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 Reference MacAyeal, Rommelaere, Huybrechts, Hulbe, Determann and RitzMacAyeal 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, a_zz – 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, H₀(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’Reference Bainess (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.

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

and

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’ = g (ρ_a −ρ)/ρ₀ 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 c_d =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.

Entrainment

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

obtained by Reference GasparGaspar (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, ΔB_a ≥ 0 is the buoyancy difference across the mixed layer/ambient ocean interface and ΔB_b ≥ 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 T_b and S_b, 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 Reference GasparGaspar, 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_* ³) /(ΔB_bṁ). 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, K_h, 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 K_h = A_h, and for the control case below we use A_h = K_h = 10m²s^-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 T₀ = -2.0°C and S near S₀ = 34.5 psu. Water at the unresolved molecular interface between ice and ocean is assumed to have properties T_b and S_b 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 (Reference Holland and JenkinsHolland and Jenkins, 2001)

Here ṁp₁c1(T_b – T₁) represents the energy needed to heat the ice from its internal temperature, T₁ , to the in situ melting point, Tb, £ = 3.35 x 10⁵ 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, Y_T and YS, are given by

where v₀ = 1.95 x 10 ~⁶ m² 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 (Reference DurranDurran, 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 (Reference Johnson, Münchow, Falkner and MellingJohnson 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. Reference BainesBaines, 2008). We cannot model this separation with the equations presented above, so we always select T_a(z) and S_a(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 Reference Holland and FelthamHolland and Feltham, 2006; Reference Holland, Feltham and JenkinsHolland 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 D₀ = 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 U² _*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 (|U₀|) represents the rms speed of the unresolved tidal flow. For all the calculations presented below, unless otherwise indicated, we use U*₀ = 0.0025 ms^-1, corresponding to an rms tidal speed of 0.05 ms^-1, given our choice of C_d. 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.

Choice of inflow ice profile, H₀(x)

A thickness profile H₀(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 https://cms.cresis.ku.edu/node/5351 in November 2011). We therefore proceed with the ‘4 plus 12’ thickness profile, as a case which exhibits a degree of realism.

Spin-up procedure

To spin-up the control case we imposed the inflow ice thickness profile, H₀(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 Reference RignotRignot 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 (Reference RignotRignot 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 (Reference RignotRignot 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, v₀), 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 H₀(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 − v₀) + v₀ 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, v₀ = 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 (Reference Johnson, Münchow, Falkner and MellingJohnson 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 m² 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 ∂²b/∂x² > 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 (∂²b/∂x² > 0).

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 km³ 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 km³ 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 km³ 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 (Reference Johnson, Münchow, Falkner and MellingJohnson and others, 2011).

Temperature perturbations

Recent warming in nearby Nares Strait (Reference Munchow, Falkner, Melling, Rabe and JohnsonMunchow and others, 2011) and the likelihood that warming of coastal waters has contributed to ice retreat and acceleration elsewhere around Greenland (Reference Holland, Thomas, de Young, Ribergaard and LyberthHolland and others, 2008; Reference MurrayMurray and others, 2010; Reference ChristoffersenChristoffersen and others, 2011; Reference StraneoStraneo and others, 2011) motivated a set of experiments in which the ambient water properties, T_a(z) and S_a(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 km³ 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.3^oC, 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.1^oC 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. 3^oC subsurface warming. Cooling the subsurface water by 0.3^oC 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 km³ 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 H₀(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 Reference GasparGaspar, 1988). We were able to obtain qualitatively similar results using the CMO scheme, which we do not discuss further here. A third scheme, from Reference Zilitinkevich and MironovZilitinkevich 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 Reference Holland and FelthamHolland 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 (Reference Johnson, Münchow, Falkner and MellingJohnson 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 m²s^-1 shuts off channel formation by smoothing and broadening the regions of intense melting. Increasing the mixed-layer minimum thickness, D₀, 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 v₀ 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 ₀ 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 5^oC 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. A_h K_h U_*0 and j) in order for the model to qualitatively match the observed geometry of Petermann Gletscher. Only A_h, and K_h, appear in other applications of the same ocean model. Reference Payne, Holland, Shepherd, Rutt, Jenkins and JoughinPayne and others (2007) and Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland and others (2009) used A_h = K_h, = 100 m² s^-1 with a grid spacing of 1 km. Reference Holland, Feltham and JenkinsHolland and others (2007) used A_h, = K_h, = 500 m²s^-1 with a grid spacing of 2 km. These are all considerably larger than our value of A_h , = K_h, = 10 m² 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.