Monte Carlo ice flow modeling projects a new stable configuration for Columbia Glacier, Alaska, c. 2020

Due to the abundance of observational datasets collected since the onset of its retreat (c. 1983), Columbia Glacier, Alaska, provides an exciting modeling target. We perform Monte Carlo simulations of the form and flow of Columbia Glacier, using a 1-D (depth-integrated) flowline model, over a wide range of parameter values and forcings. An ensemble filter is imposed following spin-up to ensure that only simulations that accurately reproduce observed preretreat glacier geometry are retained; all other simulations are discarded. The selected ensemble of simulations reasonably reproduces numerous highly transient post-retreat observed datasets. The selected ensemble mean projection suggests that Columbia Glacier will achieve a new dynamic equilibrium (i.e. “stable”) ice geometry c. 2020, at which time iceberg calving rate will have returned to approximately preretreat values. Comparison of the observed 1957 and 2007 glacier geometries with the projected 2100 glacier geometry suggests that Columbia Glacier had already discharged ∼ 82 % of its projected 1957–2100 sea level rise contribution by 2007. This case study therefore highlights the difficulties associated with the future extrapolation of observed glacier mass loss rates that are dominated by iceberg calving.


Introduction
The transfer of land-based ice into the ocean is now the leading cause of sea level rise (cf.Bindoff et al., 2007), providing almost twice the sea level rise contribution as the ther-mal expansion of sea water (∼ 55 and 30 % of total sea level rise respectively; Cazenave and Llovel, 2010).During the 1991-2002 period, small glaciers and ice caps external to the ice sheets contributed 0.77 ± 0.26 mm a −1 of sea level rise.Alaskan glaciers have the most negative total mass balance of glaciated regions outside the ice sheets (Kaser et al., 2006).A comparison of digital elevation models suggests that Alaskan glaciers contributed 0.12 ± 0.02 mm a −1 to sea level rise over the 1962-2006 period (Berthier et al., 2010).Laser altimetry observations indicate an Alaskan glacier sea level rise contribution of 0.27 ± 0.10 mm a −1 between 1992 and 2002 (Arendt et al., 2002).This latter contribution rate, however, is considered an overestimate, due to the extrapolation of glacier centerline altimetry data across glacier width.Dynamic thinning reaches a maximum along a glacier centerline, and reaches a minimum at the lateral margins of a glacier (Berthier et al., 2010).The Alaskan glacier contribution to sea level rise has also been examined in several gravimetry studies.Chen et al. (2006) inferred a contribution of 0.28 ± 0.06 mm a −1 over the 2002-2005 period.Luthcke et al. (2008) subsequently inferred a contribution of 0.23±0.01mm a over the 2003-2007period. Most recently, Jacob et al. (2012) inferred a contribution of 0.13 ± 0.02 mm a −1 over the 2003-2010 period.Together, these observations suggest that the Alaskan contribution is equivalent to ∼ 8 % of the total observed sea level rise over the 1993-2007 period (Cazenave and Llovel, 2010).
Of all Alaskan Glaciers, Columbia Glacier is presently the single largest contributor to sea level rise.Over

W. Colgan et al.: Monte Carlo ice flow modeling
the 1995-2001 period, Columbia Glacier contributed ∼ 7.1 km 3 a −1 of water to sea level rise, equivalent to ∼ 0.6 % of total observed sea level rise over the 2003-2007 period (Arendt et al., 2002;Cazenave and Llovel, 2010).Prior to the c. 1983 onset of its rapid and ongoing retreat, Columbia Glacier had an area of ∼ 1070 km 2 and a length of ∼ 66 km (Meier et al., 1985;Krimmel, 2001).The preretreat terminus position, first documented in 1794, is believed to have been stable since the fifteenth century (Rasmussen et al., 2011).Since 1983, Columbia Glacier has retreated ∼ 18 km and lost ∼ 100 km 2 of ice-covered area from its terminus (Fig. 1).This rapid retreat has been well documented, which makes Columbia Glacier an exciting modeling target (Movie 1, full movie is available in the Supplement associated with this article).Observational datasets suitable for model validation include: (i) the pre-retreat ice surface elevation profile (Meier et al., 1985), (ii) the pre-retreat ice surface velocity profile (Meier et al., 1985), (iii) the contemporary surface mass balance rate profile and mean equilibrium line altitude (Mayo, 1984;Tangborn, 1997;Rasmussen et al., 2011;O'Neel, 2012), (iv) a time-series of terminus position (Krimmel, 2001), and (v) a time-series of surface ice velocity at ζ = 50 km (Krimmel, 2001), where ζ is the curvilinear coordinate system describing downstream distance on Columbia Glacier's main flowline (complete variable notation provided in the Appendix).While not strictly an observed quantity, time-series of iceberg calving rate have also been inferred for Columbia Glacier (Krimmel, 2001;Rasmussen et al., 2011).
We examine the past and future behavior of Columbia Glacier using a 1-D (depth-integrated) flowline model that incorporates longitudinal coupling stresses and uses statistical parameterizations for two important, but poorly understood, processes: basal sliding and iceberg calving.We execute Monte Carlo simulations over a wide parameter space, to identify the cumulative uncertainty associated with both parameter and forcing uncertainties, and to provide robust ensemble mean histories and projections of variables of interest.We use an ensemble filtering technique to eliminate unrealistic simulations, whereby specific simulations are discarded if they do not: (i) satisfactorily reproduce observations of ice thickness (a state variable) at the conclusion of a transient spin-up, or (ii) initiate retreat within 100 yr of the onset of a transient forcing.Monte Carlo selection approaches have been used extensively in the context of oceanic (e.g.Van Leeuwen and Evensen, 1996) and atmospheric (e.g.Anderson and Anderson, 1999) modeling.In glaciology, Monte Carlo simulations have been used to explore uncertainty in basal sliding velocity and surface mass balance rate parameters (Chandler et al., 2006;Machguth et al., 2008;Gardner et al., 2011).
Deterministic modeling of tidewater glaciers is predicated on the implicit assumption that tidewater glaciers are intrinsically predictable and follow defined trajectories to stable attractor states, whereby small changes in initial condi-tions and/or parameters result in small changes in trajectories.Given the possibility of true chaotic behavior of tidewater glaciers, however, whereby small changes in initial conditions and/or parameters result in large changes in trajectories, the behavior observed at a given tidewater glacier is just one of a large number of possible trajectories (M.Lüthi, personal communication).By executing a large number of simulations over a wide parameter space, and then selecting simulations that reproduce observed behavior, an ensemble filtering technique provides the framework to quantify and assess non-deterministic behavior.Monte Carlo simulation also offers a powerful approach for quantifying uncertainties in observed variables resulting from uncertainties in initial conditions, parameterizations and forcing.When combined with ensemble filtering, whereby simulations are screened based on their agreement with observations, the technique can constrain the values of initial conditions, parameterizations and forcing.Although the Monte Carlo ensemble filtering is computationally intensive, it avoids the limitations of simpler uncertainty propagation approaches that often employ linearization.For the highly nonlinear problem of ice flow, which is subject to many interacting sources of uncertainty, traditional linear calculations of uncertainty propagation are unlikely to be accurate.
The stochastic probing we perform as part of the model parameter space identifies the plausible bounds of poorly constrained parameters, such as maximum accumulation rate at high elevations, while also producing thousands of simulations that plausibly explain the observed trajectory of Columbia Glacier.We find that the available diverse observational datasets are reasonably well reproduced by the ensemble mean of the selected simulations.When projected into the future, the selected ensemble mean simulation indicates that Columbia Glacier will achieve a new dynamic equilibrium geometry (i.e."stable" position), and hence no longer significantly contribute to sea level rise, by c. 2020.Thus, this case study suggests that caution must be exercised in the future extrapolation of contemporary mass loss rates that are dominated by the highly transient variations of iceberg calving rate.

Ice flow model
We apply a previously published (Colgan et al., 2012) Meier et al. (1985) to describe the "main" flowline (M) and tributaries "west" (W), "main-west" (MW) and "east" (E).Annual terminus position over the 1984 to 2010 period is also shown (updated from Krimmel, 2001).Inset: Location of Columbia Glacier in Alaska.
ice thickness (∂H / ∂t) according to mass conservation: where b is annual surface mass balance rate, w is the glacier width and ∂Q / ∂x is the along-flowline divergence of ice discharge.Following Marshall et al. (2005), depth-integrated ice discharge (Q) is taken as: where F is a spatially variable and dimensionless correction factor (discussed below), u b is the basal sliding velocity, A is the flow law parameter (we assume that Columbia Glacier is at the pressure-melting-point throughout and take A as 140 MPa 3 a −1 ; O' Neel et al., 2005), n is the flow law exponent (taken as 3), ρ is the glacier density (taken as 900 kg m −3 ), g is gravitational acceleration (taken as 9.81 m s −2 ), ∂z s / ∂x is the ice surface slope and τ is driving stress, taken as the sum of both gravitational and longitudinal coupling stresses.In an approximation of the momentum balance, depth-averaged longitudinal coupling stress (τ xx ) is included as a perturbation to the gravitational driving stress ( Van der Veen, 1987;Marshall et al., 2005): Depth-averaged longitudinal coupling stress is calculated according to Eq. ( 21) of Van der Veen (1987).This formulation derives longitudinal coupling stress by solving a cubic equation describing equilibrium forces independently at each node, based on ice geometry and prescribed basal sliding velocity.Following the suggestion of Van der Veen (1987), we assume that the longitudinal gradients of the depth-averaged longitudinal deviatoric stress are small, so that the "D" term in his Eq. ( 11) may be neglected, producing a simpler form of his Eq. ( 21), which becomes our Eq.( 4): As noted by Van der Veen (1987), this formulation is similar to the Alley and Whillans (1984) approximation for depthaveraged longitudinal coupling stress.Flowline models for alpine glaciers often invoke a parameterization to account for lateral effects on ice flow due to finite or variable glacier width (i.e.Paterson, 1994).Implementing a traditional "shape factor" parameterization of τ , however, only accounts for the influence of cross-valley shape on Q due to internal deformation, and neglects the influence of cross-valley shape on Q due to basal sliding, by implicitly assuming that basal sliding is acting on an infinitely wide glacier (i.e.Paterson, 1994).In most alpine glaciers, a shape factor parameterization of τ is valid, as internal deformation rather than basal sliding comprises the majority of Q.At Columbia Glacier, however, Q due to basal sliding is significantly greater than Q due to internal deformation throughout the ablation zone (Kamb et al., 1994;Pfeffer, 2007).Therefore, the influence of cross-valley shape on Q due to basal sliding cannot be ignored.In the spirit of a shape factor, we prescribe a spatially variable correction factor (F ), to account for the influence of cross-valley shape on both Q due to internal deformation and basal sliding.Eq. ( 2) describes ice flow within a wide rectangular cross-valley multiplied by the geometric correction factor F , and thus may be interpreted as accounting for the influence of cross-valley shape on Q due to both internal deformation and basal sliding.Incorporating w into this equation provides a rigorous expression for mass flux.We prescribe tuned, spatially variable values of F that are informed by local glacier geometry (glacier half-width divided by ice thickness; w / (2H ); Fig. 2).When calculating w / (2H ), we use the pre-retreat centerline ice thickness (H ) inferred by McNabb et al. (2012) and glacier width (w) interpolated from the distance measured between lateral shear margins along the www.the-cryosphere.net/6/1395/2012/The Cryosphere, 6, 1395-1409, 2012 main flowline of Columbia Glacier in the 1 : 100 000 Plate 5 map of Meier et al. (1985).
Making the now common assumption that the contribution of internal deformation to surface ice velocity is negligible in the ablation zone of Columbia Glacier (i.e.downstream of ∼ 40 km; Kamb et al., 1994;Pfeffer, 2007) allows us to implement a statistical parameterization of basal sliding velocity.This empirical, and hence site specific, parameterization is predicated on the observation that ice surface velocity profiles observed over the 1981 to 2001 period (Pfeffer, 2007) can be approximated with a simple exponential curve of the form: where k is a dimensional coefficient of 1 m a −1 , x is the distance downstream from km 0 and α is a scaling length (Fig. 3).This basal sliding prescription is not a sliding rule, whereby basal sliding velocity is parameterized to vary with glacier geometry or hydrology, but rather a curve fit of observed sliding velocity as a function of flowline distance (x); similar to a curve fit of surface ablation as a function of elevation (z; Eq. 6).Observations indicate that α ranged between ∼ 8.9 km in 1981 and ∼ 5.8 km in 2001, depending on terminus position.We prescribe α as a function of terminus position (x term ), which allows α to decrease as the terminus retreats upstream.The above basal sliding prescription theoretically allows basal sliding to occur anywhere along the flowline.The range of α values we impose, however, practically restrict significant basal sliding to only the ablation zone of the flowline, consistent with observations.We assume that α reaches a minimum of 5.25 ± 0.25 km when the terminus position retreats to km 50, the approximate upstream limit of the inferred bedrock over-deepening of the main flowline of Columbia Glacier (McNabb et al., 2012).The assumption that km 50 is a stable terminus position is couched in the notion that a stability criterion, comprised of the ratio between ice thickness (H ) and water depth (H w ), can distinguish stable and unstable terminus positions of tidewater glaciers.Empirical evidence suggests that tidewater terminus geometry may be regarded as stable when H / H w ≥ 1.5, and unstable when H / H w < 1.5 (Pfeffer, 2007).We use inferred bedrock elevation and observed 2007 ice surface elevation (McNabb et al., 2012) to calculate the H / H w profile along the main flowline of Columbia Glacier.These observations suggest that H / H w < 1.5 downstream of km 50, where water depth is large compared to ice thickness, but H / H w ≥ 1.5 upstream of km 50, where water depth is small compared to ice thickness (Fig. 4).Thus, we make the important assumption that the basal sliding profile will cease to evolve once the terminus retreats upstream of km 50.In each Monte Carlo simulation we randomly perturb α by a value uniformly distributed between −0.25 and +0.25 km.As α resides in an exponent, this parameter range yields a wide variety of basal sliding profiles for a given terminus position.For example, perturbing the 1992 velocity profile approximation by α = 6.8 ± 0.25 km results in an ensemble velocity range of ±1.0 km a −1 at km 55, and ±2.5 km a −1 at km 60 (Fig. 3).
Similar to Nick et al. (2007), we parameterize annual surface mass balance rate (b) as a linear function of ice surface elevation (z s ) according to: where γ is the observed annual surface mass balance rate gradient ( b / z s ; taken as 0.0085 / a; Rasmussen et al., 2011), z ela is the equilibrium line altitude, and b max is the maximum surface mass balance rate (i.e.accumulation or snowfall rate).Randomly prescribing z ela from a uniform distribution between 850 and 1050 m and b max from a uniform distribution between 3.0 and 6.0 m a −1 yields a range of surface mass balance rate profiles that encompass the empirical range (Mayo, 1984;Tangborn, 1997;Rasmussen et al., 2011;O'Neel, 2012;Fig. 5).During spin-up, z ela is prescribed as 200 m lower than the contemporary range (i.e. from a uniform distribution between 650 and 850 m) to simulate the cooler climate with which pre-retreat Columbia Glacier was most likely in equilibrium (Nick et al., 2007).

Climatic variability and forcing
In order to simulate natural climatic variability, we introduce a stochastic element by allowing equilibrium line altitude to randomly vary each decade (i.e.z ela ± δz ela ).The magnitude of the decadal perturbation (δz ela ) is randomly selected from a distribution derived from reanalysis data (Compo et al., 2011).We assume that annual z ela variability ( z ela / t) may be approximated by dividing annual air temperature variability, the difference in mean melt season air temperature from year to year ( T / t), by local environmental lapse rate ( T / z) at equilibrium line altitude: This assumes that equilibrium line altitude is correlated with a given isotherm during the melt season (e.g.Andrews and Miller, 1972).In order to determine appropriate values of T / t and T / z, we extract 137-yr time-series of 900 and 950 mb melt season (1 April to 30 September) air temperature at Columbia Glacier from Twentieth Century Reanalysis V2 Data (Compo et al., 2011;Fig. 6).The 900 mb pressure level corresponds to ∼ 990 m elevation, the approximate equilibrium line altitude of Columbia Glacier over the reanalysis period.Reanalysis data suggests that during the 1871 to 2008 period, the mean local environmental lapse rate ( T / z) was 6.7 K km −1 , and the annual variability in mean melt season air temperature ( T / t) exhibited an approximately normal distribution centered on 0 K a −1 (Fig. 6 inset).Dividing this T / t distribution by the mean local environmental lapse rate yields a distribution of annual z ela variability ( z ela / t; Eq. 7).We convert this annual z ela / t distribution into a decadal z ela / t distribution by applying a 10-yr running mean to 10 000 yr of synthetic z ela variability generated using the annual z ela / t distribution (Fig. 7).This synthetic data suggests that decadal z ela perturbations (δz ela ) can be described by a normal distribution with a mean of 0 m and a standard deviation of 30 m.
During transient spin-up, equilibrium line altitude is perturbed each decade around a fixed mean z ela .During the subsequent transient forcing period, however, the mean z ela is also forced upwards based on the long-term air temperature trend ( T / t).The long-term trend in T / t is taken as the linear trend in the 900 mb air temperature.In each Monte Carlo simulation, long-term T / t is randomly prescribed from a uniform distribution between 0.0057 and 0.0262 K a −1 .This range corresponds to the minimum and maximum trends (i.e.trend ± standard slope error) in air temperature over the 1871 to 2008 period (dashed lines Fig. 6).Dividing this rate of air temperature increase ( T / t) by local environmental lapse rate ( T / z) yields the rate of z ela increase ( z ela / t) imposed during the transient forcing period Eq. ( 7).This future climate forcing conservatively assumes no acceleration in the contemporary rate of increase in air temperature.

Model implementation and boundary conditions
We apply the 1-D depth-integrated flowline model described in Sect.2.1 (Colgan et al., 2012) to the main centerline of the Columbia Glacier.The differential equations describing transient ice thickness (∂H / ∂t) were discretized in space using first-order finite volume methods ( x = 250 m).The semidiscrete set of ordinary differential equations was then solved using ode15s, the stiff differential equation solver in MAT-LAB R2008b with a time-step ( t) of 1 yr.The numerical code does not appear to demonstrate any sensitivity to prescribed time-step over a tested range of 1/12 ≤ t ≤ 2. We selected t = 1 to facilitate the direct comparison of model output with the available observed annual datasets, without performing temporal interpolation of the model output.The model was optimized to run on eight parallel processors using the parallel computing toolbox in MATLAB R2008b.The mean processor time per Monte Carlo simulation was ∼ 48 s (Fig. 8).This allowed 20 000 simulations to be completed in ∼ 33 wall-clock hours using a 750 W Dell PowerEdge 2950 server with eight 2.83 GHz processors and a total of 32 GB of RAM.
The model ice geometry is initialized with observed preretreat ice surface elevation (Meier et al., 1985) and inferred bedrock elevation (McNabb et al., 2012).Prescribed surface mass balance rate is a source/sink term in the ice flow model.Basal sliding velocity is also prescribed externally in the ice flow model.The surface (top) boundary condition of the ice flow model, the assumption that τ → 0 at the free surface of the glacier, is implicit in the first-order formulation of the Navier-Stokes equations described by Eq. ( 3).The upstream (left) boundary condition is a second-type (prescribed flux) Neumann boundary condition to simulate an ice flow divide (i.e.Q = 0 at x = 0 km).
The downstream (right) boundary condition at the glacier terminus is a first-type (prescribed head) Dirichlet boundary condition, as the ice discharge at the terminus node (Q term ) is not known.This empirical, and hence site-specific, downstream boundary condition is based on the observation that mean terminus ice cliff height has varied between 80 and 100 m since 1981 (Pfeffer, 2007).In each simulation, the prescribed ice cliff height is randomly selected from a uniform distribution between 80 and 100 m, to assess model sensitivity.At the conclusion of a time step, terminus position is explicitly updated as the node downstream of which ice   2007(McNabb et al., 2012)).Tidewater terminus geometry may be regarded as stable when H / H w ≥ 1.5 and unstable when H / H w < 1.5 (Pfeffer, 2007).
surface elevation is less than the prescribed ice cliff height; all ice downstream of this node is prescribed to calve.While this calving parameterization honors the observed terminus ice cliff height of Columbia Glacier, we acknowledge that it is not physically based, in comparison to parameterizing calving rate as a function of longitudinal strain-rate (Nick et al., 2010).We note that an overarching goal of the Monte Carlo ensemble filter approach is to explore the response of a diverse population of Columbia Glaciers to a range of transient forcings, rather than to replicate or isolate an individual process.Thus, similar to the basal sliding and surface mass balance rate parameterizations we prescribe, a site-specific empirical calving parameterization facilitates our exploration of stable and unstable states of Columbia Glacier.Total iceberg calving rate (D) is taken as the sum of both transient ice discharge through the terminus node (Q term ) and the prescribed change in terminus position due to imposed iceberg calving: where subscript i denotes node index, and H is a Heaviside function of the form: where x crit is the location where ice surface elevation is equivalent to the prescribed ice cliff height.While the inclusion of correction factor (F ) and glacier width (w) in the calculation of ice discharge Eq. ( 2) account for flow divergence and convergence stemming from changes in glacier width, by implicitly modifying ∂Q / ∂x  Mayo, 1984;Tangborn, 1997;Rasmussen et al., 2011;O'Neel, 2012), and the parameterized range used in this study (dashed lines; Eq. 6).
with ∂F / ∂x and ∂w / ∂x terms, this parameterization does not account for the influence of tributaries.The main flowline of Columbia Glacier receives discharge from three major tributaries: "west" at ∼ km 51, "east" at ∼ km 38 and "main-west" at ∼ km 29, respectively (Fig. 1).We explicitly account for tributary effects by increasing ice inflow at the junction of each tributary by an amount proportional to the main flowline ice discharge.This additional ice inflow is smoothly distributed over several adjacent nodes using a Gaussian curve (1 km standard deviation).We increase ice inflow by temporally invariant tunable factors of 80, 25 and 40 % at km 29, 38 and 51, respectively.While these factors are imposed at tributary junctions, they represent the additional ice inflow not only from the tributary, but also the numerous smaller glaciers and cirque basins between tributaries.For example, a comparison of the pre-retreat centerline velocities of the similar-sized main and main-west tributaries (600 and 300 m a −1 , respectively; Meier et al., 1985) suggests main-west likely contributed an additional 50 % ice inflow to the main flowline at km 29.There are, however, ∼ 6 smaller glaciers/cirque basins between km 0 and 29, which we estimate to contribute the remaining 30 % additional ice discharge at km 29.

Monte Carlo ensemble filtering
We executed a large number of model simulations (20 000) in order to provide a robust ensemble mean projection of specific variables of interest, and also assess the cumulative effect of both parameter and forcing uncertainties.We randomly varied four key model parameters over a wide parameter space, two of which influence surface mass balance rate (b max and z ela ), and two of which influence ice flow (α and ice cliff height).We also randomly varied the main forcing parameter, the rate of increase in 900 mb air temperature ( T / t).Each simulation begins with a 500-yr fully transient spin-up.At the conclusion of this 500-yr spin-up, the first ensemble selection filter was imposed: only simulations that reproduced observed pre-retreat (i) mean ice surface elevation between km 40 and 60 to within ±100 m (Meier et al., 1985) and (ii) terminus position (x term ) to within ±2 km (Meier et al., 1985) were selected to carry forward into a 250yr transient forcing period.Simulations that did not satisfactorily reproduce features (i) and (ii) were discarded.The wide parameter space of the selected ensemble of simulations produced a population of modelled Columbia Glaciers of varying sensitivities (where "sensitivity" is broadly defined as mean ice reservoir overturn time in the spirit of Johannesson et al., 1989).Relatively high basal sliding and surface accumulation simulations yielded glaciers with lower mean ice reservoir overturn time than relatively low basal sliding and surface accumulation simulations.
During the subsequent 250-yr transient forcing period, this selected population of glaciers was forced by a wide range of rates of increase in equilibrium line altitude.A second ensemble selection filter was imposed to discard simulations in which retreat did not initiate within 100 yr of forcing onset.
As retreat initiated at different times between simulations, the floating model time of the twice selected simulations (i.e.those which accurately reproduced pre-retreat glacier geometry and initiated retreat within 100 yr of forcing onset), was transposed to real time by a least-squares fit between modelled and 24-yr observed terminus position histories.Subjecting the selected population of glaciers, with varying climatic sensitivities, to a wide range of climatic forcings produced a robust ensemble mean history and projection for a number of observable variables including: equilibrium line altitude, terminus position, velocity at km 50 and iceberg calving rate.The spread across the selected ensemble provides a robust measure of the cumulative uncertainty resulting from both parameter and forcing uncertainties.

Results
An inherent trade-off exists between the number of simulations selected and the size of the parameter space; a larger parameter space decreases the probability that a given simulation will achieve selection criteria but increases the robustness of the ensemble mean.Of the 20 000 Monte Carlo simulations initialized, 3022 (∼ 15 %) passed the first selection filter at the end of the 500-yr transient spin-up and were carried forward into the 250-yr transient forcing period.The remaining 16 978 simulations (∼ 85 %), which failed to reproduce observed pre-retreat ice geometry at the end of spinup, were not carried forward into the transient forcing period.Of the 3022 simulations carried forward, 353 were discarded by the second selection filter, as they did not exhibit a retreat within 100 yr of the onset of forcing.Thus, 2669 simula- tions (∼ 13 %) passed both ensemble selection filters.The selected ensemble exhibited a slight preference for terminus ice cliff height < 92 m, in comparison to ice cliff height > 93 m (Fig. 9).We regard this sensitivity as low, however, as the mean terminus ice cliff height of the selected 2669 simulations is only 2 m less than the mean ice cliff height initially prescribed to all 20 000 simulations.
The selected simulations contain the full range of initial equilibrium line altitude values (650 to 850 m) and maximum surface mass balance rate values (3.0 to 6.0 m a −1 ) over a wide range of basal sliding velocities (Fig. 10).The population of selected simulations appears to exhibit a preference for high sensitivity simulations (i.e.relatively high maximum surface mass balance rate (or accumulation rate) and basal sliding values and relatively low equilibrium line altitude) in comparison to low sensitivity simulations (i.e.relatively low maximum surface mass balance rate (or accumulation rate) and basal sliding values and high equilibrium line altitude).We note that only 5 % of the selected simulations exhibited a maximum surface mass balance rate < 4.5 m a −1 .We interpret this as the minimum high elevation accumulation rate required for sufficient mass input to maintain Columbia Glacier's pre-retreat geometry.
Both the ice surface elevation and velocity profiles of the selected simulations at the conclusion of transient spin-up, taken to be representative of the pre-retreat profiles, compare well with 1977/78 observed ice surface elevation and velocity profiles interpolated at every second kilometer along the main flowline of Columbia Glacier (Meier et al., 1985;Fig. 11).While the ensemble mean modelled velocity profile generally captures the shape of the observed velocity profile, some discrepancies exist.Firstly, the modelled profile fails to capture the localized velocity influence of an icefall at ∼ km 23.The failure of the model to adequately represent the complex physics at an icefall, where significant crevassing www.the-cryosphere.net/6/1395/2012/The Cryosphere, 6, 1395-1409, 2012 occurs, likely stems from the momentum balance approximation employed; the assumption of continuum mechanics is not valid where ice becomes discontinuous.Secondly, the modelled profile underestimates surface ice velocity in the vicinity of km 35.This is likely due to an underestimation of local convergence.This suggests that the measured distance between lateral shear zones may not be a good proxy for glacier channel width in the vicinity of km 35.Finally, the modelled ice velocity at km 66 (the terminus) slightly underestimates the velocity assessed by Meier et al. (1985).We note that the 1977/78 velocity observations downstream of ∼ km 62 are not in situ, but rather extrapolated from upstream photogrammetric values (Meier et al., 1985).
In addition to achieving good agreement with observed pre-retreat ice surface elevation and velocity profiles, the modelled ensemble mean time-series of equilibrium line altitude, terminus position, ice velocity at km 50 and calving rate also agree well with previously published observed and inferred records (Mayo, 1984;Tangborn, 1997;Krimmel, 2001;Rasmussen et al., 2011;O'Neel, 2012;Fig. 13).We note that these previously published equilibrium line altitudes represent period means, and are therefore constant over their respective time intervals, while our modelled equilibrium line altitude is transient.The combination of (i) a 200 m depression of equilibrium line altitude to simulate "cooler" climate during spin-up, and (ii) an air temperature forcing of between 0.0057 and 0.0262 K a −1 (combined with a mean local environmental lapse rate of 6.7 K km −1 ), produces an ensemble mean equilibrium line altitude that agrees well with observed contemporary equilibrium line altitude.If climate forcing persists at its current rate, the mean equilibrium line altitude at Columbia Glacier will be ∼ 1200 m by the year 2100.The ensemble mean suggests that Columbia Glacier will achieve a new stable terminus position at ∼ km 42 (near the grounding line) c. 2020, and maintain this terminus position until at least 2100 (Fig. 13).While modeled ice geometry is sensitive to the prescribed basal sliding velocity profile, we note that we prescribe a wide range of basal sliding velocity profiles to the ensemble of simulations and find the selected profiles to closely match observed profiles over the observational period.We acknowledge, however, that the future ice geometry we project is highly dependent on the assumption that the prescribed stability criterion (H / H w ) is time independent.At present, empirical evidence suggests that the stability criterion we have selected is representative of several well-studied glaciers over a wide range of time periods (Pfeffer, 2007).The modelled time-series of transient terminus position may suggest a slightly faster retreat than observed.We regard any discrepancy in retreat rate as within uncertainty across the ensemble, defined by the envelope of Monte Carlo simulations, and discuss possible causes of a slight mismatch in Sect. 4. Differences between the 1977/78 and 2100 ice surface elevation and velocity profiles are generally restricted to the region downstream of the km 35 convergence with the main-west tributary (cf.Fig. 11b, d).The absence of significant changes to ice geometry and velocity upstream of the km 23 icefall, both prior to and after retreat, is noteworthy (Fig. 12).The apparent stability in both ice geometry and velocity upstream of the icefall are direct consequences of the assumption of a time independent tidewater stability criterion.With a bedrock elevation of 1240 m above sea level at the km 23 ice fall, however, we speculate that rapid changes in basal sliding velocity associated with tidewater instability are unlikely to occur above the ice fall.
The ensemble mean time-series of surface ice velocity at km 50 generally reproduces the magnitude of velocity changes observed at km 50 (Krimmel, 2001;Fig. 13).The ensemble of simulations indicate that the surface ice velocity at km 50 increases by a factor of between 2.5 and 5.0, relative to pre-retreat (i.e. 1977/78) velocities, between the onset of retreat and the time when the terminus retreats upstream of km 50.The finer features of the km 50 velocity record, however, such as the precise timing of acceleration and temporal velocity inflections, are not well reproduced.We note that surface ice velocity essentially reflects basal sliding velocity in the ablation zone of Columbia Glacier (Kamb et al., 1994;Pfeffer, 2007), and despite employing a rather simple basal sliding parameterization, this model framework achieves a good agreement between observed and modelled ice velocity at km 50.
The ensemble mean iceberg calving rate time-series suggests that iceberg calving will "turn off" (i.e. return to dynamic equilibrium values) c. 2020, when a new stable terminus position is achieved, just as quickly as iceberg calving "turned on" at the initiation of retreat c. 1983 (Fig. 13).Thus, the total response time of Columbia Glacier to the retreat initiated by contemporary climate forcing is expected to be ∼ 40 yr.There is good agreement between ensemble mean modelled and inferred iceberg calving rate until c. 1995.After c. 1995, modelled calving rate begins to decrease, slowly until c. 2005 and then more quickly until c. 2020, while inferred calving remains elevated until at least 2007 (Rasmussen et al., 2011).This discrepancy likely stems from compounding errors during the calculation of iceberg discharge.The decrease in modelled iceberg calving rate coincides with a c. 2002 minima in both F and w (Fig. 14).Any errors in F and w are compounded when calculating iceberg calving rate Eqs.( 2) and ( 8).
The calving term also compounds uncertainty in the two statistical parameterizations used to represent basal sliding velocity and change in terminus position due to iceberg calving.While these statistical parameterizations achieve good first-order agreement with ice geometry and velocity observations, they are undeniably less robust than first-principles physically based parameterizations.Finally, part of the discrepancy between modelled and inferred calving rate is due to the fact that the inferred rate pertains to the entire Columbia Glacier complex, both the west and main branches, while the modelled calving rate only applies to the main branch once the terminus retreats upstream of the km 51 confluence.This distinction, however, should only result in discrepancy after c. 2005, when the terminus position retreats upstream of km 51.

Model limitations
While five diverse observed datasets -(i) pre-retreat ice surface elevation profile, (ii) pre-retreat ice surface velocity profile, (iii) contemporary surface mass balance rate profile and mean equilibrium line altitude, (iv) time-series of terminus position, (v) time-series of surface ice velocity at km 50 following the onset of retreat -are reasonably well reproduced by the model, the 1-D flowline model does not accurately reproduce iceberg calving rate after c. 1995.
We note that our modelled time-series are derived from a mass conserving numerical framework, unlike the diverse observed and inferred time-series.Consequently, modelled iceberg calving rate and ice velocity are interdependent.Acknowledging this interdependence, we opt to match the timeseries of observed ice velocity at km 50 at the expense of the time-series of inferred iceberg calving rate.
While the ensemble of simulations selected based on ice geometry does not appear to be sensitive to the prescribed ice cliff height, ice cliff height and calving flux are clearly related.For example, ice discharge would be 20 % greater through a 100 m ice cliff at flotation, than an 80 m ice cliff at flotation.While our modelling framework prescribes a variety of ice cliff heights across the simulations, ice cliff height is constant in a given simulation.Photographic analysis, however, suggests that terminus ice cliff height has not been constant during the retreat of Columbia Glacier (E.Welty, personal communication).The failure to acknowledge that ice cliff height reached a maximum in the Kadin-Great Nunatak (K-GN) gap at km 53 is expected to result in a proportional underestimation of calving flux during this period.The 1-D flowline model also suffers from inherent limitations in the treatment of: (i) lateral effects (i.e.convergence/divergence due to complex bed topography/tributaries) and (ii) glacier density.Complex lateral effects stemming from bed topography are a significant issue in the vicinity of the K-GN bedrock constriction at km 53.The lateral effects stemming from bedrock topography at the constriction are complicated by the lateral effects of converging ice flow from the west tributary immediately upstream (km 50 to 53).The K-GN bedrock constriction is represented in the model by a minimum glacier width of 3 km prescribed at km 53, based on the distance between lateral shear margins in the 1977/78 ice surface velocity map (Meier et al., 1985;Fig. 14).The aerial photography record has revealed the emerging bedrock topography of the K-GN gap exposed by thinning ice.The total pre-retreat glacier width of 5 km at km 53 has now declined by 60 %, to just 2 km between the bedrock shores of the now ice-free K-GN gap (Fig. 1).Changes in glacier width over the retreat period are not as pronounced elsewhere along the flowline.While employing a transient correction factor, i.e.F (w(t),H (t)) rather than a constant correction factor, i.e.F (w,H ) may offer some potential to refine the treatment of a bedrock constriction in a flowline model, it would not improve the treatment of tributary convergence.A 2-D (plan view) model offers a better potential to improve the treatment of the bedrock constriction than further parameterization of a 1-D (flowline) model.Generally, however, even with 1-D limitations of lateral effects, the ice geometry and timing of retreat is reasonably well reproduced as the glacier retreats through the bedrock constriction.
Similar to previous Columbia Glacier modeling investigations (O'Neel et al., 2005;Nick et al., 2007) we assume that glacier density is constant in space and time (taken as 900 kg m −3 ).At Columbia Glacier, however, observations suggest that heavy crevassing can result in extremely low bulk glacier densities in the ablation zone (e.g.< 700 kg m −3 in the top 85 m of ice at km 63.7; Meier et al., 1994).Furthermore, these observations, as well as anecdotal evidence, suggest that the ablation zone of Columbia Glacier has become progressively more crevassed since the retreat began c. 1983 (Meier et al., 1994).Continuity calculations suggest that glacier density decreases by ∼ 20 % as ice flows downstream from the bedrock constriction at km 53 to the glacier terminus, achieving depth-averaged bulk glacier densities as low as 750 kg m −3 (Venteris, 1997).Thus, in reality, glacier density at Columbia Glacier is neither constant in time nor space.This has important consequences for an ice flow model predicated on mass conversation with invariant density.For example, an increase in bulk glacier density over time would result in an increase in ice volume over time, which would decrease the apparent modelled rate of terminus retreat (i.e. the upstream migration of the terminus due to calving would be offset by the volumetric expansion of remaining ice).Temporally variable glacier density, stemming from changes in crevasse spacing or extent, is expected to influence both surface mass balance rate and ice dynamics (Colgan et al., 2011).Spatially and temporally transient glacier density, however, is not incorporated in even the most sophisticated ice flow models, including Elmer (Gagliardini and Zwinger, 2008), Community Ice Sheet Model (CISM; Lipscomb et al., 2009) and Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009).It is not immediately apparent how to derive an appropriate equation of state, or even statistical parameterization, that would allow rate of change of glacier density to be incorporated into the mass continuity equation.
Time-lapse photography of Columbia Glacier's flow just upstream of the K-GN gap at km 53 provides a compelling visualization of the complex flow we are modelling with a 1-D flowline model (Movie 1).The time-lapse photography, compiled from Extreme Ice Survey photographs obtained over the 12 May 2007 to 18 April 2012 period, especially highlights the prevalence of crevasses in the Columbia Glacier ablation zone.Given the evident complexity of ice flow just upstream of the K-GN gap, it is encouraging that the 1-D flowline model framework reasonably reproduces the observed terminus retreat rate through the reach documented by Movie 1, as well as the ice velocity observed at km 50, just upstream from the camera position.The time-lapse photography also illustrates the complexity of the discrete iceberg calving events being parameterized.We find that time lapse photography provides not only unique insight to the entirety of complex glaciological processes, but also qualitative validation of model parameterizations employed to capture the form and flow of Columbia Glacier during its highly transient retreat.

Projecting sea level rise
Forecasting the cryospheric contribution to sea level rise over the next century is a task of paramount importance for the glaciology community.At present, climatically-forced projections of the small glacier and ice cap contribution to global sea level rise have been restricted to the surface mass balance component (Raper and Braithwaite, 2006;Radiae and Hock, 2011), ignoring the potential contribution due to ice dynamics (primarily due to the paucity of data required to initialize and validate ice flow models, and secondarily due to the computational expense associated with solving the transient equations for ice flow for a large population of glaciers).Rule-based projections of land ice mass change, including, but not limited to, statistical extrapolations, have been applied to both ice sheets (e.g.Velicogna, 2009;Rignot et al., 2011) and small glaciers and ice caps (e.g.Meier et al., 2007;Pfeffer et al., 2008).These projections extrapolate future sea level rise contributions from trends contained in the observed record, and include both climatically-controlled surface mass balance changes and dynamic-controlled changes in calving rate.Our present study of Columbia Glacier serves as a reminder that caution should be exercised when performing such projections (e.g.Price et al., 2011).While the absolute magnitude of our modelled iceberg calving rate does not agree precisely with observations after c. 1995, the good agreement between observed and modelled ice geometry and velocity lends high confidence to the general timing and shape of iceberg calving rate projection.Thus, it is very likely that the iceberg calving rate of Columbia Glacier will indeed "turn off" (i.e. return to a stable dynamic equilibrium value) c. 2020 just as quickly as it "turned on" c. 1983 (Figs. 12 and 13).This exemplifies how a statistical extrapolation of pre-2000 trends into the post-2000 period could lead to erroneous projections.As iceberg calving rate is not a "state" variable (i.e. a variable capable of predicting the future behavior of a system), extrapolating past mass loss rates that are primarily based on iceberg calving is not robust, as rapid changes in ice dynamics can result in decade-scale transitions between stable and unstable states.
We estimate the remaining response of Columbia Glacier to the rapid retreat initiated c. 1983 by comparing the total sea level rise contributions of the model domain in the 2007-2100 and 1957-2100 periods.We can calculate the total anticipated sea level rise contribution from the model domain (i.e. the main flowline of Columbia Glacier) between 2007 and 2100 by quantifying the difference in ice volume (δV ) between the inferred 2007 ice thickness (H 2007 ;McNabb et al., 2012) and the ensemble mean modelled 2100 ice thickness (H 2100 ), according to This formulation projects a total sea level rise contribution of 16.0 km 3 from the main flowline model domain over the 2007 to 2100 period.This estimate excludes all ice-covered areas outside the main flowline model domain (i.e.adjacent small glaciers, cirques and tributaries).only a fraction of the total ice-covered area of the Columbia Glacier complex (∼ 260 of ∼ 910 km 2 ; McNabb et al., 2012), these numbers are gross underestimates of the absolute values of the sea level contribution of the entire ice complex (which was 160 km 3 over the 1957-2007 period;McNabb et al., 2012).They do, however, illustrate that the majority of the response of Columbia Glacier to the terminus perturbation initiated c. 1983 has been completed.While this exercise is similar to the concept of "committed sea level rise" (Price et al., 2011), it differs slightly by maintaining a climate forcing throughout the entire transient simulation.

Summary remarks
We apply a 1-D (depth-integrated) flowline model to Columbia Glacier that incorporates longitudinal coupling stresses and statistical parameterizations for basal sliding and  Mayo, 1984;Tangborn, 1997;Krimmel, 2001;Rasmussen et al., 2011;O'Neel, 2012) time-series of equilibrium line altitude (z ela ; A), terminus position (x term ; B), ice surface velocity at km 50 (u 50 s ; C) and inferred calving flux (D; D) at Columbia Glacier over the 1850 to 2100 period.iceberg calving.A computationally efficient implementation allows Monte Carlo simulations to be executed over a wide parameter space to produce robust histories and projections of variables of interest, as well as assess the cumulative effect of both parameter and forcing uncertainty.Ensemble selection filters are imposed at: (i) the conclusion of spin-up, to ensure an accurate reproduction of pre-retreat glacier geometry, and (ii) 100 yr into the forcing period, to ensure terminus retreat has initiated.The resultant twice selected ensemble of simulations reproduces several observed datasets within the uncertainty envelope defined by the ensemble range.Inferred iceberg calving rate is not well reproduced in the vicinity of the K-GN bedrock constriction at km 53.A 2-D (plan-view) model is required to resolve the complexities of ice flow in this region.It is not clear, however, how the issue of significant glacier density transience, if indeed occurring, may be resolved with a modeling approach predicated on a continuum mechanics momentum balance (i.e. the assumption that ice does not become discontinuous at any time or place).
The ensemble mean projection suggests that Columbia Glacier will achieve a new stable ice geometry c. 2020, by which time iceberg calving rate will have decreased to a dynamic equilibrium value much lower than that observed during the highly transient 1990s and 2000s.Comparison of the pre-retreat (1957) and 2007 glacier geometries with the projected 2100 glacier geometry suggests that, by 2007, Columbia Glacier had already discharged ∼ 82 % of the total sea level rise contribution expected by 2100.As the model  As Columbia Glacier is the single biggest contributor to sea level rise from Alaska, knowledge of its individual response to highly transient contemporary climate change is useful for predicting the response of the Alaskan tidewater glacier population.Based on the sheer magnitude of pre-retreat Columbia Glacier, in terms of ice-covered area (∼ 1070 km 2 ), ice thickness (up to 975 m) and bedrock overdeepening (up to 525 m), it has been suggested that Columbia Glacier may be considered analogous to a marine-based ice sheet (Molnia, 2008).Thus, understanding and predicting the response of Columbia Glacier to contemporary climate change likely also has value in anticipating the response of the ice sheets.
depthintegrated (1-D) flowline model with a first order approximation for longitudinal coupling stresses to the main centerline of Columbia Glacier.The model domain of the center flowline of Columbia Glacier extends from the main flow divide at ∼ 2750 m elevation at km 0 (61.369• N and 147.153 • W) down to sea level at km 70 (60.974• N and 147.093 • W; Fig. 1).The model solves for the transient rate of change in

Fig. 2 .
Fig. 2. Observed pre-retreat ratio of glacier half-width to ice thickness (w / (2H )) along the centerline of Columbia Glacier, with the corresponding spatially variable correction factor (F ) applied to the ice flow model in this study.

Fig. 3 .
Fig. 3. Observed ice surface velocity (u s ) profiles at Columbia Glacier over the 1981 to 2001 period (solid lines; Pfeffer, 2007) and their corresponding parameterizations (dashed lines; Eq. 5) using differing values of exponential length scale (α).Grey shading denotes α ± 0.25 km around the 1992 profile.Inset: The empirical relation between exponential sliding length scale (α) and terminus position (x term ) used in this study.

Fig. 4 .
Fig. 4. The ratio between observed ice thickness (H ) and water depth (H w ) along the main flowline of Columbia Glacier in2007(McNabb et al., 2012)).Tidewater terminus geometry may be regarded as stable when H / H w ≥ 1.5 and unstable when H / H w < 1.5(Pfeffer, 2007).

Fig. 6 .
Fig. 6.Mean melt season (1 April to 30 September) 900 mb air temperature (T ) over the 1871 to 2008 period at Columbia Glacier extracted from the Twentieth Century Reanalysis V2 data provided by NOAA/OAR/ESRL PSD(Compo et al., 2011).Inset: Corresponding histogram and non-parametric distribution of annual variability in 900 mb air temperature ( T / t).

Fig. 7 .
Fig. 7. Synthetic annual ( t = 1 a) variability in equilibrium line altitude ( z ela ) over 10 000 yr, generated using the T / t distribution shown in Fig. 6 and a lapse rate ( T / z) of 6.7 K km −1 .The corresponding decadal and centurial variability are also shown ( t = 10 and 100 a, respectively).Inset: Histogram and normal distribution (mean = 0 m; standard deviation = 30 m) of decadal z ela perturbations (δz ela ).

Fig. 8 .
Fig. 8. Histogram of processor time per simulation of the 20 000 Monte Carlo simulations.The dashed line denotes the ensemble mean (48 s).The bimodal distribution is due to the greater computational requirements of simulations selected to carry forward into transient forcing following spin-up (B) in comparison to those that were not selected (i.e.discarded following spin-up; A).

Fig. 9 .
Fig. 9. Prescribed terminus ice cliff height in the selected ensemble of simulations.Dashed line denotes the ensemble mean (88 m).

Fig. 10 .
Fig. 10.Equilibrium line altitude (z ela ), maximum surface mass balance rate (b max ; or accumulation rate) and mean basal sliding velocity between km 50 and 60 (u b ) in the selected ensemble of 2669 simulations.Dashed line denotes the fifth percentile of accumulation rate of the selected ensemble of simulations.

Fig. 11 .
Fig. 11.Modelled (grey lines with ensemble mean in black) and observed (points; Meier et al., 1985) ice surface elevation (z s ) and ice surface velocity (u s ) along the Columbia Glacier main flowline (x) in 1977/78 (pre-retreat; A and B) and in 2100 (post-retreat; C and D).

Fig. 12 .
Fig. 12. Modelled time-space evolution of ensemble mean rate of change of ice thickness (∂H / ∂t) and surface ice velocity (u s ) along the Columbia Glacier main flowline (x) between 1970 and 2100.Colorbars saturate at −100 m a −1 and 2000 m a −1 , respectively.The black line denotes the ensemble mean terminus position over the period.

Fig. 14 .
Fig. 14.Ensemble mean time-series of terminus glacier width (w term ) and correction factor (F term ).
suggests a short response time (∼ 40 yr) between the initial perturbation of the tidewater terminus (c.1983) and the attainment of a new dynamic equilibrium geometry (c.2020), this case study highlights the difficulties associated with extrapolating glacier mass loss estimates that are dominated by iceberg calving into the future.