High sensitivity of tidewater outlet glacier dynamics to shape

Variability in tidewater outlet glacier behavior under similar external forcing has been attributed to differences in outlet shape (i.e., bed elevation and width), but this dependence has not been investigated in detail. Here we use a numerical ice flow model to show that the dynamics of tidewater outlet glaciers under external forcing are highly sensitive to width and bed topography. Our sensitivity tests indicate that for glaciers with similar discharge, the trunks of wider glaciers and those grounded over deeper basal depressions tend to be closer to flotation, so that less dynamically induced thinning results in rapid, unstable retreat following a perturbation. The lag time between the onset of the perturbation and unstable retreat varies with outlet shape, which may help explain intra-regional variability in tidewater outlet glacier behavior. Further, because the perturbation response is dependent on the thickness relative to flotation, varying the bed topography within the range of observational uncertainty can result in either stable or unstable retreat due to the same perturbation. Thus, extreme care must be taken when interpreting the future behavior of actual glacier systems using numerical ice flow models that are not accompanied by comprehensive sensitivity analyses.


Introduction
While recent dynamic changes in marine-terminating outlet glaciers in Greenland and Antarctica are broadly correlated to climatic and oceanographic conditions, substantial spatiotemporal variability is evident (e.g., Howat et al., 2008;Moon and Joughin, 2008;McFadden et al., 2011;Moon et al., 2012;Walsh et al., 2012).Glaciers in close proximity, and presumably exposed to similar environmental forcing, display contrasting behavior, suggesting that their dynamic response is largely dependent on individual characteristics, such as glacier shape (Meier and Post, 1987;Howat et al., 2007Howat et al., , 2008;;Pfeffer, 2007).
The dynamics of tidewater glaciers are sensitive to changes in the stresses acting on their marine termini that result from changes in external forcing.Any reduction in resistance to flow due to mélange weakening (Amundson et al., 2010;Howat et al., 2010;Christoffersen et al., 2011), ice tongue thinning (Holland et al., 2008;Motyka et al., 2011) and/or grounding line retreat (Joughin and Alley, 2011) leads to acceleration and a transfer of resistive stresses to the glacier margin and bed.This acceleration results in increased stretching and thinning (i.e., dynamic thinning).For glaciers with relatively shallow slopes that are grounded well above flotation, thinning will reduce gravitational driving stress and discharge, causing the ice to slow and stabilize.If, however, the ice thins to near the point of flotation, thinning will likely reduce friction at the ice bed, causing additional acceleration (Pfeffer, 2007).This positive feedback can lead to unstable, runaway retreat if the glacier is grounded across a basal depression, as retreat of the grounding line into deeper water will further increase the grounding line discharge, resulting in rapid thinning to flotation within the depression (Schoof, 2007).As shown in Gudmundsson et al. (2012), an increase in lateral ice flow convergence can limit this positive feedback and stabilize the grounding line prior to unstable retreat across a basal depression.However, outlet glaciers that lack lateral ice flow convergence along their topographically confined trunks should be particularly sensitive to the aforementioned positive feedback between retreat, acceleration, and E. M. Enderlin et al.: High sensitivity of tidewater outlet glacier dynamics to shape thinning across a basal depression, as demonstrated by the rapid retreat of Helheim Glacier through a basal depression from 2004-2005(Howat et al., 2007;;Nick et al., 2009).
The timing and total magnitude of retreat will therefore depend on the basal topography and changes in glacier width, as rises in the bed and lateral constrictions in the surrounding bedrock walls should act as stabilizing points of ice convergence and higher friction (O'Neel et al., 2005;Jamieson et al., 2012).The response of a glacier to a perturbation at its front, therefore, should be highly dependent on the shape of the valley through which it flows (Pfeffer, 2007;Jamieson et al., 2012).Understanding this dependence is important for assessing regional variability in glacier behavior, identifying glaciers likely to exhibit large-scale changes in the near future, and constraining the impact of measurement uncertainty on model predictions of glacier behavior.

Model description
To test the influence of valley shape on tidewater glacier dynamics, we perform sensitivity tests using a movinggrid, depth-integrated, width-averaged numerical ice flow model (Vieli and Payne, 2005;Nick et al., 2009Nick et al., , 2010;;Vieli and Nick, 2011) that includes lateral, basal, and along-flow stresses and uses an effective pressure-dependent sliding law and crevasse depth-dependent calving law (Benn et al., 2007;Nick et al., 2010).The depth integration of the model implicitly employs the Shallow Shelf Approximation, which is not fully appropriate for the entire model domain.However, the model results obtained from the sensitivity tests examined herein should be valid along the regions of parallel flow within the topographically confined trunks of fastflowing tidewater outlet glaciers, as demonstrated by similar type model's ability to reproduce changes in observed tidewater glacier behavior in numerous studies (see Nick et al., 2009Nick et al., , 2012;;Vieli and Nick, 2011;Colgan et al., 2012).Details on the shape and surface mass balance parameterizations, governing equations, boundary conditions, and the applied perturbation are provided below.The model discretization and implementation procedures are described in detail in Appendix A.

Shape and surface mass balance parameterizations
Each model glacier consists of a 120 km-wide inland accumulation basin (ice sheet) that drains into a narrow, topographically confined outlet channel with a bed below sea level.Within the outlet, we assess the effects of both width and width gradient along-flow using six different width configurations that are within the range commonly observed for fast-flowing tidewater outlet glaciers in Greenland (Fig. 1a) and an idealized, over-deepened bed profile based on available bed elevation measurements (Fig. 1b).The cross-sectional shape of the outlet channel is also likely to vary along-flow and could be incorporated into the model through the use of a shape factor, but its influence on ice flow is outside of the scope of the width-averaged sensitivity tests examined herein.
The bed elevation profiles are based on measurements for Helheim Glacier, Kangerdlugssuaq Glacier, and Jakobshavn Isbrae in Greenland (CReSIS, http://www.cresis.ku.edu/data).We applied best-fit polynomials to along-flow bed elevation profiles from each glacier in order to extract elevation and slope ranges, which were then used to construct an idealized bed (Fig. 1b, solid black line).A basal depression was included in the idealized bed profile because (1) similar depressions were found in all three bed elevation maps and (2) glaciers overlying basal depressions should be particularly sensitive to force balance perturbations, as described above.Using the uniform 7 km-width profile, we then assess the effects of bed elevation uncertainty using three additional bed profiles that fall within the ≤ 50 m-uncertainty of current ice thickness observations (acquired from radio echo-sounding) (Bamber et al., 2013) as shown in Fig. 1b: shoal depth decreased by 35 m (dashed-dotted line), depression depth decreased by 35 m (dashed line), and shoal and depression depths decreased by 35 m (dotted line).
The surface mass balance (SMB) rate is held constant in time and is prescribed as a function of distance from the equilibrium line.The magnitude of accumulation varies slightly for the different outlet shapes in order to maintain a similar interior ice thickness, as would be observed for glaciers fed by the same catchment area.The resulting SMB profiles fall within the typical range for Greenland outlet glaciers (Ettema et al., 2009;Burgess et al., 2010).
We include submarine melting along the base of the floating ice when present.Submarine melting is temporally constant but varies spatially as a function of distance from the grounding line (Rignot and Steffen, 2008), with the maximum melt rate of 0.6 m d −1 occurring ∼ 1.2 km from the grounding line.Submarine melt rate magnitudes are based on the range of melt rates estimated for west Greenland outlet glaciers in Enderlin and Howat (2013).

Governing equations
Assuming no ice flow transverse to the glacier flowline, the temporal change in ice thickness can be determined using conservation of mass, such that where H is ice thickness, U is the vertically averaged horizontal ice velocity, W is width, B is surface mass balance (including submarine melting), t is time, and x is distance from the ice divide along the central flowline.The assumption of parallel flow is valid for our purposes, since we are focusing on dynamic changes near the front where ice flow is confined by shear margins and/or fjord walls.The governing force balance equation determined through conservation of momentum is where β 2 is the basal friction coefficient, A is the rate factor, ρ i = 917 kg m −3 is the density of ice, g is gravitational acceleration, h is the ice surface elevation, and v is the depthaveraged effective viscosity, which is defined as The RHS of Eq. ( 2) is the gravitational driving stress, which is balanced by gradients in longitudinal stress (1st term LHS), basal resistance (2nd term LHS), and lateral resistance (3rd term LHS).The rate factor, A, is scaled with the cumulative strain rate (A k ∝ k j =1 ∂U ∂x j ) in order to account for strain heating along-flow.Using this scaling, the rate factor increases from a minimum of 3.5 × 10 −25 Pa −3 s −1 at the divide to a maximum of 1.7 × 10 −24 Pa −3 s −1 at the calving front, corresponding to a depth-averaged ice temperature range of −10 • C to −2 • C (Cuffey and Paterson, 2010).The basal friction coefficient, β 2 , is parameterized as the product of the basal roughness factor and basal effective pressure.The same decreasing, piecewise linear function is used to parameterize the basal roughness factor for all simulations.
We assume an open connection between the ocean and icebed interface, such that the water pressure increases with the bed depth and the basal effective pressure equals zero at the grounding line.The values selected for β 2 are similar to those used to model Helheim Glacier (Nick et al., 2009), with values decreasing from a maximum of β 2 ≈ 7.0 × 10 10 Pa s m −1 to zero as the ice approaches flotation.
The grounding line location is tracked using a flotation criterion, which has been successfully used in similar models to reproduce observed grounding line migration (e.g., Nick et al., 2009;Jamieson et al., 2012).The model employs a moving-grid that adjusts the grid spacing, x, at each time step to precisely and continuously track the location of the grounding line by stretching/contracting the coordinate system to maintain x ∼ = 200 m (see Appendix A for details).

Boundary conditions
The up-glacier boundary is the ice divide (i.e., U = 0) and the down-glacier boundary (i.e., calving front) is located at the point along the ungrounded portion of the glacier trunk where the surface crevasse depth equals the ice surface elevation (Benn et al., 2007;Nick et al., 2010).The crevasse depth (d crev ) is calculated as follows: where ρ w = 1000 kg m −3 is the density of fresh water, d w = 10 m is the crevasse water height, and R xx is the resistive stress, which is defined as follows: At the calving front, the gradients in longitudinal stress are in balance with the difference in hydrostatic pressure between the ice and sea water such that where ρ sw = 1028 kg m −3 is the density of sea water.

Model simulations
The simulated glaciers are initialized from the same speed and thickness profiles during a 200 yr spin-up period.The speed and thickness profiles obtained at the end of the spinup are used as the initial, steady-state profiles for the perturbation experiments.Here, we define steady-state based on inter-annual variability in the ice thickness at each grid cell: a simulated glacier is considered to have reached steady state if the ice thickness at each grid cell varies by < 0.1 m yr −1 .All simulations meet this criterion by the end of the 200 yr initialization period.The transient simulations (i.e., perturbation experiments) are initialized using their respective steady-state speed and elevation profiles (Fig. 1c and d).A step reduction in resistive stress is applied at the calving front in order to simulate a reduction in backpressure resulting from ice tongue thinning and breakup, grounding line retreat, or mélange weakening.This perturbation is applied at model time step k by increasing horizontal stretching (i.e., decreasing resistance to horizontal flow) at the front by a factor, S, equivalent to where is the difference in hydrostatic pressure between the ice and sea water.The stress perturbation, , is the same for all simulations.By defining S in terms of the stress perturbation, we can express the perturbation in terms of an equivalent volume of ice tongue retreat (i.e., reduction in non-hydrostatic backstress).In our experiments, a constant value of = 1.00 × 10 8 Pa m is applied for the entire 30 yr-duration of the transient simulations.The magnitude of the perturbation is equivalent to the loss of up to ∼ 20 km 3 of floating ice from the terminus, which is of similar magnitude of the recent disintegration of Jakobshavn Isbrae's floating ice tongue (Joughin et al., 2004).

Model results
We focus our analysis of model results on the first 15 yr, following the application of the step perturbation when the magnitude of the simulated glacier response is largest (see Figs. 2-4).Application of the step perturbation results in instantaneous acceleration and retreat for all model runs (Figs. 2 and 3), reaching a maximum rate of thinning ranging from 11-17 m yr −1 near the grounding line and decreasing to ∼ 1 m yr −1 ∼ 35-55 km inland.Thinning and acceleration cause the discharge through the grounding line to peak within 6 months, increasing ∼ 5 % for the glaciers with narrower profiles and ∼ 10 % for the two glaciers with the widest profiles, then gradually stabilize (Fig. 4c).Following this initial response, the evolution is bimodal: for the narrower and narrowing-inland glaciers, thinning and acceleration decline from their initial increase towards a new steady state with little overall retreat (Fig. 4a and b) and increase in ice discharge (Fig. 4c), whereas the two glaciers with the widest outlets reach flotation above the basal depression, triggering a much larger retreat and discharge increase.
The initial thickness profile determines the mode of response to the perturbation, as the initial steady-state thickness profiles of wider glaciers are closer to flotation above the basal depression and have a shallower surface slope in order to maintain the same discharge across the grounding line as narrower glaciers (Fig. 1c).Less initial thinning is therefore required to reach the threshold for unstable retreat, resulting in the ungrounding of a large section of the trunk.The delay in the onset of unstable retreat (i.e., lag time) is also controlled by the initial ice thickness above the basal depression, which varies by ∼ 15 m here due to differences in convergence for glaciers that widen inland relative to those with parallel sides (Fig. 1c).Although the glacier that widens inland is initially thicker and therefore requires more time to thin to flotation (i.e., lag time), once unstable retreat is triggered, the total retreat and increase in discharge is of greater magnitude due to the feedback between retreat and increased cross-section of flow (Fig. 4).
Raising the shoal and depression elevations by 35 m in the uniform 7 km width scenario (Fig. 3, bottom right) results in less than 5 m of difference between initial surface elevations, yet it drastically influences the glacier's dynamic response.Raising the bed within the depression (Fig. 3, bottom left) increases the ice surface above flotation and reduces thinning rates by up to ∼ 5 m yr −1 in the depression, so that unstable retreat does not occur.Raising the elevation of the shoal without raising the depression (Fig. 3 thickening within the depression, quadrupling the lag time between the perturbation and unstable retreat.Once unstable retreat is triggered, however, the increase in discharge and retreat of the grounding line into the basal depression is the same magnitude as for the glacier with the original, deeper shoal (Fig. 4).Thus, small differences in the thickness gradient appear to have a similar, but stronger effect, as differences in the width gradient in controlling the timing of retreat.

Discussion and implications
All simulated glaciers that undergo unstable retreat into the depression show a significant lag time, of at least several years, between the perturbation and the onset of unstable retreat, determined by the time needed to thin to flotation above the basal depression.This result is consistent with observations from Columbia Glacier, Alaska, where a similar lag time between the onset of thinning and retreat following a period of glacier stability was observed in the 1980s (O'Neel et al., 2005).In southeast Greenland, however, a one-or twoyear lag time between elevated ocean surface temperatures and the onset of rapid, unstable retreat has been observed (Howat et al., 2008).Differences in simulated and observed lag times are likely a consequence of starting the simulated glacier from an initial steady state, whereas glaciers in south- east Greenland may have been thinning since the mid 1990s (Krabill et al., 1999;Rignot et al., 2004).Simulated lag times are also likely to be influenced by the use of a simple flotation criterion to track grounding line migration rather than the more physically based contact problem described in Nowicki and Wingham (2007), although the flotation criterion has been successfully used in similar models to reproduce observed grounding line migration (e.g., Nick et al., 2009;Jamieson et al., 2012).Further, our model shows that small variations in width and basal topography can impart large differences in the timing of unstable retreat, which may explain intra-regional variability found in Greenland.These effects can be non-local, with inland variations in width and bed elevations influencing the stability of the grounding line on decadal timescales.
A notable result of this study is that variations in ice thickness and basal topography of a magnitude similar to the accuracy of airborne radar sounding observations predict vastly different behaviors for topographically confined glaciers that are initially near flotation (i.e., within a few tens of meters) over basal depressions.It is therefore unclear whether observational capabilities are adequate for constraining prognostic simulations of such glaciers.Compounded on this problem, the spatial resolution (1-2. of the interpolated bed elevation map available for Greenland (Bamber et al., 2013) may result in smoothing of the bed topography in coastal regions, potentially influencing the simulated response (Durand et al., 2011).Our results suggest that a similar problem may exist for width, as ∼ 1 km of variation in outlet width may cause stable or unstable response.Thus, simulations of topographically confined outlet glaciers with termini near flotation must be accompanied by comprehensive sensitivity analyses to establish confidence in predictions.Further, we suggest that similar sensitivity analyses should be completed using two-or three-dimensional models in order to assess the influence of glacier shape on grounding line stability for glaciers and ice streams with strong lateral convergence along their trunks.

Conclusions
Using a simple ice flow model applied to archetypal outlet shapes, we have confirmed that the dynamic response of glaciers under a given perturbation at the ice front is highly sensitive to along-flow variations in shape, shedding light on the high spatial and temporal variability observed in outlet glacier behavior.The response of a glacier overlying a basal depression is bimodal; a perturbation results in either a grad-ual return to a new steady state with little thinning and retreat or triggers run-away, multi-kilometer retreat and tens to hundreds of meters of thinning.Whether or not a glacier will enter unstable retreat is dependent on its minimum thickness above flotation at the onset of the perturbation, which is in turn dependent on shape.For glaciers draining the same interior catchment, glaciers with wider steady-state grounding lines and those with deeper basal depressions will tend to be closer to flotation in the depression than narrower or shallow glaciers, and thus less dynamic thinning will be required to bring the ice within the depression to flotation.
Our sensitivity tests also suggest that for glaciers that are initially near flotation (i.e., within a few tens of meters) over basal depressions, both the ice thickness and bed elevation within the depression must be precisely known in order to be able to reasonably determine future stability or instability in response to external forcing.This raises the question as to whether or not current observational capabilities are adequate for constraining prognostic models of glacier behavior, as small errors may lead to substantially different predictions for glaciers that are near flotation.Although the one-dimensional numerical ice flow model utilized herein relies on several approximations that are not required in the more complex three-dimensional numerical models, the major findings of our sensitivity study are numerically robust and governed by the physics of ice flow.Thus, based on these sensitivity tests, we conclude that extreme care must be taken when analyzing numerical model results applied to actual glacier systems.

Model discretization and implementation procedures
The general discretization of the model is described in detail below.The complete Matlab ® version of the model used in this paper can be obtained by contacting the corresponding author.
Several parameters must first be specified for use throughout the model, including ice density (ρ i = 917 kg m −3 ), ocean water density (ρ sw = 1028 kg m −3 ), and gravitational acceleration (g = 9.8 m s −2 ).The initial grid spacing ( x 0 ) is used to construct the gridded length (x = 0 : x 0 : length(x)) of the model domain.The choice of x 0 should be based on desired model resolution and computation time, and was selected in this study as 200 m.At each x, temporally fixed values for bed elevation and width (h bed and W , respectively) are prescribed.Estimates for the ice thickness (H ) and surface velocity (U ) must also be specified at each grid cell for model initialization.
Equation ( 2) is used to determine the new gridded velocities for the model domain through iterative convergence.The partial differential and discretized forms of Eq. ( 2) are as follows: where v is effective viscosity defined as follows: and subscripts are position indices.Equation (A1) describes the force balance between gradients in longitudinal stress (1st term LHS), basal drag (2nd term LHS), lateral drag (3rd term LHS), and gravitational driving stress (RHS).For proper convergence to occur, the longitudinal stress term must be calculated on the staggered grid, as indicated by the position indices in Eq. (A1).The lateral drag term must be linearized so that the equation can be written in matrix-vector form.The linearization procedure is where γ j = U −2/3 j for simplification.The matrix-vector form of Eq. (A1) becomes where and the subscript c denotes the calving front location and the subscript term denotes the end of the ice domain.The calving front is located at the first ungrounded grid cell where surface crevasses generated by longitudinal stretching intersect sea level (Eqs.4-6) and the end of the ice domain is located where the ice thickness reaches zero.The boundary conditions have already been incorporated in Eq. (A4) so that there is zero ice flux at the ice divide (U 1 = 0) and the gradients in longitudinal stress are in balance with the difference between the hydrostatic pressure of the ice and ocean water at the calving face, such that The calving front boundary condition is applied from calving face to the end of the ice domain in order to avoid the force imbalance that occurs for h/ x = ∞.Ice-free grid cells are not included in the matrix-vector notation because their driving stress and velocity terms are equal to zero.
The new velocities are calculated by taking the inverse of the sparse tridiagonal coefficient matrix, M, in Eq. (A4) and multiplying by the RHS matrix (U = M −1 T ).Alternatively, the use of the Matlab backslash operator can be used to solve for U (U = M/T ) with decreased computation time and increased numerical stability.The new velocities are used to recalculate the velocity gradient and effective viscosity at each www.the-cryosphere.net/7/1007/2013/The Cryosphere, 7, 1007-1015, 2013 grid cell.This process is repeated iteratively until the difference between the velocities calculated in consecutive iterations meets a prescribed tolerance.
The gridded velocities are used to determine the change in ice thickness using conservation of mass (Eq.1), such that: where the subscript t denotes the time and t = 0.001 yr = 31 536 s is the time step.Using the results from Eq. (A7) and surface mass balance (including submarine melting), B, the new ice thickness at each grid cell is solved using: H j,t+1 = H j,t + H j,t + B j,t t. (A8) For a fixed-grid numerical model, the gridded velocities obtained for time t can be used with the ice thickness values for time t +1 (Eq.A8) to determine the coefficient matrix for time t+1.Moving-grid numerical models must account for the change in the grounding line position between time t and time t + 1 in order to accurately model grounding line migration (Pattyn et al., 2012).For the moving grid, the interp1 function in Matlab ® can be used to interpolate ice thickness values between grid cells, allowing grid adjustment to the precise location where H meets the flotation criterion.To maintain an adjusted grid spacing similar to a target value, x 0 , the new grid spacing, x t+1 , can be calculated as follows: x gz,t+1 round x gz,t+1 / x 0 , (A9) where x gz,t+1 is the location of the grounding line and "round" specifies that the divisor is rounded to the nearest integer value.The interp1 function is then used to interpolate all variables that vary spatially (e.g., H, h, U, A, etc.) to the new grid spacing.The use of relatively small grid spacing ensures that no numerical diffusion is introduced into the model during this interpolation.Further, the continuous moving-grid grounding line tracking described in Eq. ( A9) is in-line with the boundary theory of Schoof (2007) and can be used to accurately capture transient grounding line migration (Pattyn et al., 2012).The interpolated variables are then used as input for Eqs.(A1-A8) to solve for gridded ice thickness and velocity at time t + 2.

Fig. 1 .
Fig. 1.Profiles of (a) width, W , and (b) bed elevation, h bed , along the outlet channel.Steady-state profiles of (c) surface elevation, h, and (d) speed, U , obtained for each glacier shape.The different width and bed profiles used to obtain the steady-state profiles in (c) and (d) are distinguished by the line colors and styles specified in (a) and (b), respectively.

Fig. 2 .
Fig.2.Modeled annual elevation and speed profiles within the outlet channel (color coded for time, see colorbar) obtained using the 6 width configurations examined herein for 15 yr following the onset of the step perturbation.In the elevation profiles, the thin dashed gray line is the ice surface elevation required to remain grounded.

Fig. 3 .
Fig.3.Modeled annual elevation and speed profiles within the outlet channel (color coded for time, see colorbar) obtained using the 4 bed elevation profiles examined herein for 15 yr following the onset of the step perturbation.In the elevation profiles, the thin dashed gray line is the ice surface elevation required to remain grounded.As in Fig.1, the different bed elevation profiles are distinguished by line style.

Fig. 4 .
Fig. 4. Time series of modeled (a) front position change, (b) grounding line position change, both relative to their respective initial positions, and (c) the fractional increase in grounding line discharge following the perturbation onset relative to the initial steady-state discharge.