The sensitivity of flowline models of tidewater glaciers to parameter uncertainty

. Depth-integrated (1-D) ﬂowline models have been widely used to simulate fast-ﬂowing tidewater glaciers and predict change because the continuous grounding line tracking, high horizontal resolution, and physically based calving criterion that are essential to realistic modeling of tidewater glaciers can easily be incorporated into the models while maintaining high computational efﬁciency. As with all models, the values for parameters describing ice rheology and basal friction must be assumed and/or tuned based on observations. For prognostic studies, these parameters are typically tuned so that the glacier matches observed thickness and speeds at an initial state, to which a perturbation is applied. While it is well know that ice ﬂow models are sensitive to these parameters, the sensitivity of tidewater glacier models has not been systematically investigated. Here we investigate the sensitivity of such ﬂowline models of outlet glacier dynamics to uncertainty in three key parameters that inﬂuence a glacier’s resistive stress components. We


Introduction
Width-and depth-integrated (1-D) numerical ice flow models (i.e., flowline models) have been used to simulate the dynamic behavior of narrow, fast-flowing tidewater glaciers in a number of geographic settings, including Greenland (e.g., Nick et al., 2009Nick et al., , 2012Nick et al., , 2013Vieli and Nick, 2011), Antarctica (e.g., Gladstone et al., 2012;Jamieson et al., 2012), Iceland (Nick et al., 2007a), and Alaska (Nick et al., 2007b;Colgan et al., 2012). These models simulate ice flow by balancing the gravitational driving stress with basal and lateral resistive stresses and along-flow longitudinal stress gradients. The parameterization of basal and lateral resistive stresses vary for each model: basal resistance is commonly described by a Weertman-type sliding law that assumes an effective-pressure dependency (e.g., Nick et al., 2009;Vieli and Nick, 2011;Jamieson et al., 2012) and lateral resistance is parameterized by integrating the horizontal shear stress over the channel width (van der Veen and Whillans, 1996) or by multiplying the driving stress by a shape factor that accounts for differences in the cross-sectional area of the glacier along flow (Cuffey and Paterson, 2010, pp. 342).
The inclusion of longitudinal stresses accounts for the along-flow transfer of stress perturbations, such as those arising from ice shelf thinning and grounding line retreat. In these cases, accelerated ice flow near the terminus will be transferred inland through longitudinal stress coupling, leading to dynamic thinning along the outlet glacier trunk (Nick et al., 2009;Vieli and Nick, 2011). In order to accurately simulate both the initial change in dynamics triggered by external forcing, and the resulting evolution, the grounding line position must be accurately tracked using moving grid or heuristic approaches that satisfy the Schoof (2007) boundary layer theory . Changes in the calving front position on seasonal (or longer) timescales must also be determined using a physically based calving criterion (e.g., Nick et al., 2010;Vieli and Nick, 2011) rather than an arbitrary, fixed calving front position. Although several large-scale ice sheet models include longitudinal stress gradients and continuous grounding line tracking Gudmundsson et al., 2012;Cornford et al., 2013), most large-scale models are currently unable to incorporate physically based calving front variations necessary for simulating rapid changes in tidewater glacier dynamics.
All models must be constrained by observations. Although ice thickness and speed data are available for many tidewater glaciers, little or no data are available on the properties of the ice (i.e., temperature, fabric, damage, etc.) and interactions between the ice and the underlying bed, leading to large uncertainty in the appropriate parameter values for the ice rheology and the basal sliding relation. These parameters are, therefore, either tuned until the simulated ice thickness and speed reasonably approximate the available measurements (e.g., Nick et al., 2009Nick et al., , 2013 or solved using inverse methods (e.g., Morlighem et al., 2010;Arthern and Gudmundsson, 2010), usually under the assumption of initial stability. While variations in these parameters within their uncertainty have been shown to strongly influence ice sheet model predictions (Stone et al., 2010;Greve et al., 2011;Applegate et al., 2012), the sensitivity of tidewater glacier models to these parameters remains largely untested, with most studies focusing on uncertainty in the external forcing (Vieli and Nick, 2011;Colgan et al., 2012;Nick et al., 2013).
Here we examine the sensitivity of a moving-grid, widthand depth-integrated numerical ice flow model (Nick et al., 2009(Nick et al., , 2010(Nick et al., , 2013Vieli and Nick, 2011) to uncertainty in three key ice rheology and basal sliding parameterizations. While recent modeling has shown that, in some cases, transverse flow can strongly effect the behavior of tidewater glaciers , the simplicity of flowline models allow for more straightforward assessment of the relative sensitivity to individual parameters and their efficiency allows for exploration of a large area of parameter space. We also note that for a typical tidewater glacier flowing through a rock-walled fjord and terminating at a grounded front or confined ice tongue, the assumption of parallel flow inherent in flowline models is likely to be valid.
In Sect. 2, we describe the geometry of the simulated glaciers, the key parameters included in our sensitivity study and their influence on ice flow, as well as the external forcing parameterizations used in our model simulations. Section 3 describes the results of the initial stable and transient (i.e., perturbed) simulations. Namely, we find that a nonunique combination of parameter values can produce similar (i.e., within the typical measurement uncertainty) stable glacier configurations that exhibit dramatically different responses to a step change in external forcing. The influence of parameter uncertainty on the model results is presented in Sect. 4 and the implications for flowline models applied to real glacier systems are presented in Sect. 5. We find that varying model parameters within observational constraints can result in vastly different predictions of glacier behavior from similar initial conditions and conclude that, in potentially many cases, prognostic models of tidewater glaciers cannot be reasonably constrained using the available data.

Model description
We use a previously published ) flowline model that includes lateral, basal, and along-flow longitudinal stresses and uses an effective pressure-dependent sliding law as well as a crevasse depth-dependent calving criterion (Benn et al., 2007;Nick et al., 2010) to test flowline model sensitivity to parameter uncertainty. The simulated ice flow depends on a number of parameters, including ice rheology, basal sliding, and surface mass balance. Here we focus on two parameters used to describe the ice rheology, namely the rate factor and enhancement factor, as well as the exponent used to relate ice flow speed and basal resistive stress (i.e., basal sliding relation exponent). We focus our analysis on these three parameters in particular because they strongly influence the resistive stress balance components that govern ice flow but are poorly constrained for most fast-flowing tidewater outlet glaciers. The parameters are varied across 18 simulations, as described in Table 1.

Glacier geometry
The response of a tidewater glacier to external forcing is strongly influenced by the basal topography through the feedback between ice thickness and the discharge (i.e., the volume of ice passing across the grounding line per unit time, Pfeffer, 2007;Schoof, 2007). Although an increase in lateral ice flow convergence can limit this positive feedback and stabilize the grounding line on a reverse bed slope for ice streams , this stabilizing mechanism is absent for outlet glaciers confined by bedrock walls along their lateral margins (i.e., topographically confined), making them susceptible to unstable retreat across reverse The Cryosphere, 7, 1579-1590, 2013 www.the-cryosphere.net/7/1579/2013/ The maximum width-and depth-integrated ice temperature (T max ) used to calculate the rate factor is also shown.

Model parameters and associated uncertainties
The influence of the rate factor, enhancement factor, and basal sliding relation exponent on ice flow and their prescribed uncertainty ranges are described below.

Rate factor
The rate factor, A, affects the speed at which ice deforms (i.e., strain rate) under a given stress, with higher values corresponding to lower effective viscosities and faster deformation rates. The deformation rate strongly increases with temperature, T , as dislocations in the ice become more mobile. The temperature dependence of the rate factor can be described by a simple Arrhenius relationship, with values increasing by a factor of 5-10 between 10 C and the melting temperature (Cuffey and Paterson, 2010, pp. 64).
To account for variations in ice flow due to uncertainty in the temperature-dependent rate factor, we initialize the model simulations with three climate-based temperature profiles. The temperature is initially set to a maximum of either 2, 4, or 6 C at the grounding line, decreasing inland by ⇠ 5 C per kilometer increase in surface elevation to account for colder temperatures in the ice sheet interior. During the 200 yr model spin-up, we allow the temperature to freely adjust in response to deformational heating from lateral shearing (width-integrated), frictional heating from contact between the ice and surrounding bedrock walls (widthand depth-integrated), and cooling from advection of colder ice from the interior. At each time step, the ice temperature is used to obtain the rate factor.
For the 4 C range in prescribed temperature, the rate factor at the grounding line varies from ⇠ 7 ⇥ 10 25 to ⇠ 2 ⇥ 10 24 Pa 3 s 1 , increasing non-linearly with temperature (Table 1). The temperature range spans the 5 C temperature assumed for Helheim Glacier and Jakobshavn Isbrae using similar type models (Nick et al., 2009(Nick et al., , 2013Vieli and Nick, 2011) and is extended to include higher values in order to account for potential warming of the ice in response to future changes in external forcing.

Enhancement Factor
The enhancement factor, E, is a non-dimensional scalar to the rate factor that is used to account for additional ice deformation that cannot be explained by the rate factor alone, likely due to the presence of impurities or anisotropic fabric development within the ice (Cuffey and Paterson, 2010, pp. 71). The rate factor and enhancement factor influence the effective viscosity, v, or "softness" of the ice in the same way through v = (EA) 1/3 @U @x where @U/@x is the longitudinal strain rate. Thus, larger values for either A or E result in less viscous, softer ice and www.the-cryosphere.net/7/1579/2013/ The Cryosphere, 7, 1579-1590, 2013 faster flow rates. While the rate and enhancement factors are similar in their effect, we treat them separately here because they are used to parameterize independent properties of the glacier ice (i.e., temperature and anisotropy, respectively). Further, although the rate factor can typically be constrained by measured or modeled ice temperatures, the anisotropy of the ice is often unknown. As such, in simulations of real glacier systems, the enhancement factor is tuned independently in order to reproduce measured strain rates (Cuffey and Paterson, 2010, p. 71).
The values for the enhancement factor are likely to vary by at least an order of magnitude throughout the glacier (Cuffey and Paterson, 2010, Table 3.5, p. 77), however, we only consider depth-integrated values of E = 1 and E = 2 in our model (Table 1). While laboratory experiments show that much larger enhancement factors are possible, this range is consistent with those typically used in depth-integrated flow models (Nick et al., 2009(Nick et al., , 2010(Nick et al., , 2013Vieli and Nick, 2011), and is meant to reflect the mean enhancement over the ice column.

Basal Sliding Exponent
We assume that the basal shear stress, ⌧ b , is a function of along-flow variations in the difference between the ice overburden and water pressures at the ice-bed interface (i.e., effective pressure), N e , the basal roughness factor, , and the depth-integrated ice velocity as described by where m is the basal sliding exponent (van der Veen and Whillans, 1996; Vieli and Payne, 2005;Nick et al., 2009Nick et al., , 2010. The water pressure and basal roughness factors are held constant throughout the model simulations. We assume an open connection between the ocean and ice-bed interface, such that the water pressure increases with the bed depth and N e = 0 at the grounding line. Although our effective pressure parameterization does not account for seasonal changes in subglacial water pressure, which have been shown to influence ice flow over short (i.e., hourly to monthly) scales (Bartholomew et al., 2010;Howat et al., 2010), similar type models have successfully simulated long-term changes in glacier behavior using the same parameterization (e.g., Nick et al., 2009;Vieli and Nick, 2011). The basal roughness factor is parameterized as a piecewise linear function that decreases with distance from the ice divide (Figs. 1, 2). The values for the basal roughness factor are tuned for each combination of parameter values in order to minimize the difference in ice thickness and calving front position between the simulated glaciers. We choose sliding exponent values of m = 1, 2, and 3 for our model simulations (Table 1), which are consistent with theoretical values of the friction law (e.g., Gudmundsson, 1997;Schoof, 2005) and with values used in previously published flowline studies of real glacier systems (e.g., Nick et al., 2009;Vieli and Nick, 2011;Jamieson et al., 2012).

Initial stable and transient forcing
Surface mass balance (SMB) is parameterized as a piecewise linear function of surface elevation relative to the equilibrium line altitude (ELA), which is held constant in time (Figs. 1, 2). A similar parameterization was utilized by Nick et al. (2007) and Colgan et al. (2012) for Columbia Glacier, Alaska. For the different parameter combinations, the increase in the accumulation rate with elevation above the ELA varies by ±19 % relative to the mean accumulation profile in order to maintain a similar interior ice thickness for all model simulations. In the ablation zone, the melt rate profile is the same for all model simulations, with maximum melt rates of ⇠ 1.8 m yr 1 where the glacier surface approaches sea level. The SMB profile for each simulation falls within the typical range for Greenland outlet glaciers (Ettema et al., 2009;Burgess et al., 2010). Additionally, the range of accumulation rates applied to the models is in-line with the 17 % uncertainty in SMB from the Regional Atmospheric Climate Model v.2 for the Greenland Ice Sheet (RACMO2/GR) (Ettema et al., 2009) and is likely to fall within historical variations in SMB for most glacier systems.
We model submarine melting along the base of the floating tongue as a function of distance from the grounding line (Rignot and Steffen, 2008), with the maximum melt rate of 0.4 m d 1 occurring ⇠ 1.2 km from the grounding line. During the model spin-up the submarine melt rate is held constant.
The transient model behavior is initiated by instantaneously doubling the submarine melt rate along the entire submerged ice face. Retreat of the grounding line can also be initiated, however, by maintaining the same melt rate but by shifting the maximum melt rates closer to the grounding line (Gagliardini et al., 2010). Here, we maintain the elevated rate of submarine melting throughout the transient simulation (i.e., step perturbation) to examine differences in glacier behavior to a sustained change in external forcing. Both the initial stable and transient melt rates fall within the range of estimated melt rates for Greenland floating ice tongues .

Model results
In this section, we present results for the initial stable (Sect. 3.1) and transient (Sect. 3.2) model simulations.

Stable simulations
We find that by tuning the SMB in the accumulation zone and the basal roughness factor along the entire glacier length, similar stable glacier configurations can be simulated for all 18 parameter combinations (Figs. 3, 4) the difference in ice thickness between the individual simulations and the ensemble mean (not shown) is within the up to 50 m-uncertainty of current ice thickness observations acquired from radio-echo sounding (Bamber et al., 2013). The calving front position varies by  1 km relative to the ensemble mean and flow speed at the ice front varies by less than 25 %, both of which fall within the range of observed seasonal variability in Greenland Schild and Hamilton, 2013). Thus, we consider all stable profiles to fall within the uncertainty range of the ensemble mean.

Transient simulations
The step-increase in submarine melt rate causes the ice shelf to immediately thin and the front to consequently retreat, reducing the resistive stress near the grounding line. The reduction in resistive stress leads to acceleration and thinning within the outlet channel (Figs. 5, 6). We find that the magnitude of the dynamic response varies for each simulation as the result of differences in the bed geometry and the parameter values; therefore, to isolate the influence of parameter uncertainty on the transient glacier response, the simulations performed using the different bed geometries are presented separately.

Seaward-dipping bed profile simulations
The response of the simulated glaciers to the step increase in the submarine melt rate is strongly controlled by the value of the enhancement factor: doubling the enhancement factor from E =1 to E =2 causes the change in grounding line position and discharge to approximately double following the onset of the perturbation. The pattern of retreat and acceleration following the onset of the step increase in submarine melting is the same for the E = 1 and E = 2 simulations, therefore, only the results obtained for the enhanced ice deformation simulations are presented in detail. For these simulations, we find that the response to increased submarine melting is relatively small but varies with the choice of parameter values. Within the first year of the transient simulations, the magnitude of front retreat varies by a factor of 1.5 (2-3 km) (Figs. 5, 7). The initial dynamic acceleration caused by the retreat of the floating tongue increases the volume of ice passing across the grounding line per unit time (i.e., discharge), which varies by a factor of 3 (5-15 %) for the different simulations (Fig. 7). Dynamic thinning within the outlet channel leads to grounding line retreat into shallower water by the end of the transient simulations (Fig. 7), causing the discharge to gradually decrease and adjust towards pre-perturbed values. The magnitude of grounding line retreat also varies by a factor of 3 between simulations, from a minimum of 2 km to a maximum of 6 km. The small, intra-annual variations in the extent and discharge of each simulated glacier shown in Fig. 7 are due to the use of discrete time steps in the model simulations, which cause the grounding line position to fluctuate by 10s of meters. These variations do not influence the overall response of the glacier to the perturbation.

Over-deepened bed profile simulations
In the absence of enhanced ice deformation (i.e., E =1), the simulated glaciers remain grounded seaward of the basal depression despite measurable thinning and acceleration within the outlet channel (not shown). Thinning and acceleration cause the grounding line to retreat 1.5-2.4 km into shallower water towards the crest of the shoal. Similar to the E = 1 simulations performed using the seaward-dipping bed pro-file, the grounding line discharge initially increases by 1.5-5 %, then slowly decreases towards pre-perturbed values. For the simulations with enhanced ice deformation (i.e., E = 2), the transient response to the increase in the submarine melt rate is bimodal: the glaciers either initially thin and accelerate but remain grounded on the shoal and approach a new stable configuration, or the initial thinning and acceleration bring the ice to flotation in the basal depression, triggering much larger retreat, acceleration, and thinning (Fig. 6). The magnitude of grounding line retreat and increase in discharge vary by a factor of ⇠ 1.5 (1.7-2.7 km) and 5 (0.8-4 %), respectively, for the glaciers that remain grounded across the basal depression (Fig. 8). In contrast, for the four glaciers that unstably retreat ⇠ 22 km across the basal depression, the fractional increase in discharge varies by a factor of ⇠ 2.2 (22-48 %) (Fig. 8). The

Discussion
The results of our model simulations suggest that a nonunique combination of parameter values can produce similar stable glacier configurations, as previously demonstrated by Arthern and Gudmundsson (2010). If the accumulation rate uncertainty is restricted to 5 %, the range of parameter combinations used in the model simulations can be reduced to 8 of the 18 parameter combinations. Although the culling procedure is likely to eliminate incorrect parameter combinations, the dynamic response of the culled simulations varies with the choice of parameter values used to describe the ice rheology and basal sliding.
We find that dynamic thinning rates within the outlet channel vary with the rate factor or enhancement factor of the simulated glacier, leading to differences in grounding line retreat following the onset of the step perturbation. The rate of dynamic thinning increases as the rate and enhancement factors increase, such that the magnitude of grounding line retreat is larger for glaciers with higher ice temperatures (Figs. 7,8;temperatures distinguished by color) or enhanced deformation. As demonstrated by the simulations performed using the over-deepened bed profile, differences in the ice flow from ±2 C uncertainty in the maximum ice temperature can strongly influence the dynamic response of glaciers that are initially close to flotation across reverse bed slopes because a small perturbation can bring the ice to flotation (i.e., thickness threshold), triggering a much larger response (Figs. 6, 8). For the simulations obtained with a linear basal sliding relation and enhanced ice deformation (Figs. 6,8;solid lines), dynamic thinning brings the two simulated glaciers with higher ice temperatures (T max = 2 C and 4 C) to flotation within the depression < 5 yr after the onset of the perturbation, triggering the rapid and unstable retreat of the grounding line. In contrast, the simulated glacier with the lowest ice temperature profile (T max = 6 C) remains grounded on the shoal throughout the 20 yr transient simulation (Figs. 6,8;solid blue lines). Given that the respective differences in the initial ice thickness (< 30 m) and speed (< 0.1 m d 1 ) for the T max = 4 C and T max = 6 C simulations are within the typical observational uncertainty for real glacier systems, either combination of parameter values could be used to simulate the initial stable glacier configuration in the absence of precise temperature measurements.
Differences in the basal sliding parameterization have a relatively small (< 1 m d 1 ) influence on the stable glacier configuration within the outlet channel because basal resistive stress provides little resistance to ice flow where the basal water pressure is high (Eq. 2). The choice of the basal sliding exponent does, however, strongly influence the magnitude of ice flow acceleration following the onset of the perturbation. For our model simulations, changes in basal resistive stress are proportional to changes in U 1/m , meaning that in order to increase the basal resistive stress by the same magnitude, ice flow acceleration must be larger for simulations with larger values for the basal sliding exponent. The influence of the basal sliding exponent on the transient glacier behavior is clearly shown in the discharge time series for the E =2 and T max = 2 C simulations performed using the over-deepened bed geometry ( Fig. 8; red lines). Further, the choice of the basal roughness parameter influences both the timing and magnitude of retreat, such that the response of the simulated glaciers does not necessarily vary systematically with the ice rheology, as demonstrated by differences in the grounding line retreat time series in Figs. 7 and 8.

Implications for modeling of real glacier systems
Although our simulations were performed using a depthintegrated flowline model for analytical simplicity and computational efficiency, the ice rheology and basal sliding parameterizations we test are inherent in all numerical ice flow models. Therefore, our results suggest that in some cases, prognostic ice flow models cannot be used to confidently predict the precise timing and magnitude of changes in glacier behavior due to uncertainty in the ice rheology and basal sliding parameterizations. Confidence in prognostic models can be improved, however, by comparing the modeled transient behavior with precise measurements of the initial and transient glacier configurations (i.e., hindcasting, Aschwanden et al., 2013). Using a hindcasting approach, simulations that fail to reproduce the observed temporal evolution of the glacier are rejected, restricting the range of parameters used in prog-nostic simulations. The benefit of hindcasting is clearly illustrated in Fig. 8: if similar transient results were obtained for a real glacier system with annual front position and speed time series, these data could be used to assess the validity of the different model simulations, which would likely improve the accuracy of prognostic simulations.
Of particular importance is that small differences in the rheology of the ice from temperature uncertainty of ±2 C can determine whether or not a glacier initially grounded across a basal depression will reach the threshold for unstable retreat. The high sensitivity of predicted behavior to differences in ice rheology has two major implications for modeling of real glacier systems. First, given that the uncertainty in the width-and depth-integrated ice temperature for real glacier systems is probably much larger than the ±2 C considered here, predictions of glacier behavior should be accompanied by sensitivity analyses that utilize a range of ice temperatures. Second, if the positive feedback between accelerated ice flow, heating within the lateral and basal shear margins, and ice softening within the shear margins is also considered, the temperature-sensitivity demonstrated in our model simulations calls into question the reliability of widthintegrated models. As shown with a similar flowline model applied to Jakobshavn Isbrae in western Greenland, lateral variations in ice rheology can be incorporated into a widthintegrated model by applying an enhancement factor to the lateral resistive stress component (Vieli and Nick, 2011). In that study, the inclusion of the enhancement factor improved The Cryosphere, 7, 1579-1590, 2013 www.the-cryosphere.net/7/1579/2013/ the model's ability to simulate the period of rapid grounding line retreat and acceleration. The model, however, was unable to reproduce the continued acceleration of the glacier as the rate of grounding line retreat decreased (Vieli and Nick, 2011), potentially indicating that the aforementioned feedback between acceleration and shear softening could not be accounted for using a constant enhancement factor. Thus, for real glacier systems, it is likely that temporal changes in the ice rheology and basal sliding parameterizations must also be accounted for in order to accurately simulate temporal changes in glacier behavior.
In order to improve confidence in predictions of glacier behavior, we suggest that the initial ice rheology is constrained by temperature and surface strain rate measurements that can be used to calculate the effective viscosity of the ice. Further, if available, past changes in dynamics should be used to constrain the basal sliding exponent and temporal changes in the ice rheology due to strain heating and damage. To restrict the choice of appropriate parameter values, model simulations should be performed using a range of parameter values and simulations that successfully reproduce the initial and transient glacier thickness and speed profiles, as well as the position of the calving front. Simulations that successfully reproduce past glacier configurations can then be used to predict the most likely range of glacier behavior in response to changes in external forcing.

Conclusions
In order to investigate the sensitivity of prognostic flowline models to uncertainty in ice rheology and basal sliding parameterizations, we vary the rate factor, enhancement factor, and basal sliding exponent across 18 simulations performed using a previously published flowline model. We find that although a non-unique combination of parameter values can be used to produce similar initial stable configurations, differences in the ice rheology, from uncertainty in the rate factor (through temperature) or the enhancement factor, strongly influence the magnitude of grounding line retreat in response to external forcing. Similarly, differences in the basal sliding exponent control the magnitude of glacier acceleration required to restore grounding line stability following the onset of retreat, with larger basal sliding exponents corresponding to larger fractional changes in ice flow speed. Precise measurements of the initial and transient glacier configurations can be used to restrict the range of parameter values used in the model simulations through hindcasting, however, a nonunique combination of ice rheology and basal sliding parameterizations may still exist. Based on these results, as well as the high sensitivity of flowline models to uncertainty in geometry parameterizations , we conclude that in order to confidently predict glacier behavior, flowline models applied to real glacier systems must be accompanied by sensitivity tests that take into account the un-certainty in model parameters. In the absence of such sensitivity tests, we suggest that the use of prognostic flowline models is restricted to the prediction of long-term trends rather than to precisely constraining the timing of future changes in dynamics.