Environmental controls on the thermal structure of alpine glaciers

Water entrapped in glacier accumulation zones represents a significant latent heat contribution to the development of thermal structure. It also provides a direct link between glacier environments and thermal regimes. We apply a two-dimensional mechanically-coupled model of heat flow to synthetic glacier geometries in order to explore the environmental controls on flowband thermal structure. We use this model to test the sensitivity of thermal structure to physical and environmental variables and to explore glacier thermal response to environmental changes. In different conditions consistent with a warming climate, mean glacier temperature and the volume of temperate ice may either increase or decrease, depending on the competing effects of elevated meltwater production, reduced accumulation zone extent and thinning firn. For two model reference states that exhibit commonly-observed thermal structures, the fraction of temperate ice is shown to decline with warming air temperatures. Mass balance and aquifer sensitivities play an important role in determining how the englacial thermal regimes of alpine glaciers will adjust in the future.


Introduction
Glacier ice can be cold or temperate, as defined relative to the pressure melting point.Numerous studies employing borehole thermometry (e.g.Paterson, 1971;Blatter and Kappenberger, 1988) and ice-penetrating radar surveys (e.g.Holmlund and Eriksson, 1989;Gusmeroli et al., 2010) have documented the thermal structure of glaciers.Observed thermal regimes span a range from entirely cold to entirely temperate, with different polythermal structures in between (see Irvine-Fynn et al., 2011).Theoretical and numerical studies focused on understanding the controls on and evolution of alpine glacier polythermal structure (e.g.Blatter and Hutter, 1991;Aschwanden and Blatter, 2005) have been few relative to those addressing the thermal structure of ice sheets (e.g.Robin, 1955;Dahl-Jensen, 1989;Greve, 1997b;Breuer et al., 2006;Aschwanden et al., 2012).Thermal structure is relevant to glacier hydrology (Wohlleben et al., 2009;Irvine-Fynn et al., 2011), rheology (Duval, 1977) and mass balance (e.g.Delcourt et al., 2008).An understanding of thermal structure in smaller ice masses is important for predicting their responses to changing environmental conditions (e.g.Radić and Hock, 2011).
The impacts of previous glacier states on the thermal structures of Arctic glaciers have been explored using numerical methods by Delcourt et al. (2008) and Wohlleben et al. (2009).In some cases (e.g.Rippin et al., 2011), thermal disequilibrium has been proposed as an explanation for observed thermal structure.These results lead to questions about the nature of transient thermal states and how thermal structure will evolve in the future.
Models of glacier flow often neglect the presence of water in temperate ice, despite the distinct rheological and hydrological implications (Greve, 1997a).We refer to these as temperature-based models.True polythermal models account for latent heat storage, for example, by tracking water content and freezing fronts.When applied to the Greenland Ice Sheet, Greve (1997b) finds that a polythermal model predicts a thinner layer of ice at the pressure melting point than a temperature-dependent model.Accounting for the water content of temperate ice, therefore, changes the simulated thermal structure.
Furthermore, it has not been clearly established how thermal structure evolves in a changing climate.In the future, will some glaciers become colder as suggested by Rippin et al. (2011), or will the cold ice regions in these polythermal glaciers shrink as appears to be occurring elsewhere (e.g.Pettersson et al., 2007;Gusmeroli et al., 2012)?Some Published by Copernicus Publications on behalf of the European Geosciences Union.controls on thermal structure such as surface temperature and accumulation area extent may change quickly, over years or decades.Basal heat fluxes and the rate of strain heating may adjust more slowly.Geothermal fluxes are not likely to be affected.Predictions of future glacier behaviour depend on the contributions made by multiple heat sources which may be affected by climate in dissimilar ways.
Our goal is to develop a better understanding of how changes in environmental conditions and ice dynamics shape the thermal structure of mountain glaciers.We use models applied to synthetic glacier profiles to isolate the direct influences of environmental variables in the absence of irregular or complex geometry.In the following, glacier thermal structure refers to the spatial distribution of englacial heat.Where ice is cold, this influences temperature, and where ice is temperate (at the melting point), this influences water content.The addition and removal of heat may also change the distribution of cold and temperate ice.The specific objectives of this study are (1) to evaluate the relative contributions of individual heat sources to glacier-wide thermal structure, (2) to relate the sensitivity of steady-state thermal regimes to internal and environmental variables and interpret these sensitivities with respect to observed thermal structure in alpine glaciers, and (3) to simulate the transient evolution of glacier thermal structure in response to prescribed changes in climate.

Modelling approach
We use a simple two-dimensional mechanically-coupled thermal model to both calculate steady states and to evolve thermal structure forward in time.An alternative way of representing polythermal conditions in glacier models is the enthalpy gradient method, proposed by Aschwanden and Blatter (2009).Representing thermal evolution within ice using only a single state variable (enthalpy) rather than both temperature and water content simplifies energy conservation and the model representation: separate grids for temperature and water content and explicit jump conditions between cold and temperate ice domains are not required (Aschwanden and Blatter, 2009).The capability to model a wide variety of thermal structures is also more simply implemented.The theory behind the model has been outlined in depth by Aschwanden et al. (2012).We briefly describe the model below and then present details of its implementation.

Heat flow
The flow of heat q within cold ice can be described by Fourier's law, where k is thermal conductivity and ∇T denotes a temperature gradient.Following Aschwanden et al. (2012), we replace the temperature gradient by a material enthalpy gradient ∇H (J kg −1 ) given by with specific heat capacity c p .Energy conservation results in the advection-diffusion equation where ρ is material density (either firn, ice or a combination), u is material velocity and Q is a heating rate.All terms in Eq. ( 3) have SI units of J m −3 s −1 .The diffusivity, κ, is different for cold ice (κ c ) and temperate ice (κ t ) and depends on density ρ.We parameterise the thermal diffusivity in the cold porous near-surface layer based on results by Sturm et al. (1997) as in units of W m −1 K −1 .Snow and firn are better thermal insulators than ice.
In temperate ice, temperature gradients arise only from the small pressure dependence of the melting point, so diffusive heat transfer (the first term on the right-hand side of Eq. 3) becomes negligible (Aschwanden et al., 2012).This vanishing term can be represented by setting κ t to zero, however, we follow Aschwanden et al. (2012) in prescribing κ t as a small positive constant as a means of regularisation.We choose a value two-orders of magnitude lower than κ c , small enough that the numerical solution is insensitive to further decreases in κ t .We represent the transition between κ c and κ t as a smooth function over a small enthalpy range H trans .This is done to improve numerical consistency, but also crudely represents a finite boundary layer of the type suggested by Nye (1991).

Heat sources
The source term Q in Eq. ( 3) is modelled as the sum of an internal heat source and a surface heat source: where Q str is strain heating and Q m is the heat associated with meltwater entrapment and possible refreezing.An additional heat flux (Q b ) is present at the ice-bed interface, and is the sum of the geothermal flux (Q geo ) and the frictional heat flux from basal sliding and water flow.The geothermal flux (Q geo ) component of Q b is poorly constrained in many mountainous regions as well as below the major ice sheets.We take Q geo = 55 mW m −2 as a reference value broadly representative of continental heat flux (Blackwell and Richards, 2004) and set the minimum value of Q b to this value of Q geo .The maximum value tested for Q b (1000 mW m −2 ) is larger than recent estimates of maximum continental heat fluxes by a factor of about five (Davies and Davies, 2010).Heat derived from frictional heating or dissipation from subglacial drainage has been inferred to increase the basal heating term by a factor of ten, to roughly 540 mW m −2 (Clarke et al., 1984).For the purposes of this study, we prescribe Q b directly and assume it is constant in space and time.We later justify this choice with results from sensitivity tests indicating that basal heating has a limited influence on temperate ice volume.We do not consider basal ablation because we expect it to be relatively small in most settings compared to surface ablation (cf.Alexander et al., 2011).Following Cuffey and Paterson (2010, Ch. 9), the strain heating term within a unit volume is with deviatoric stress τ and strain rate ˙ .Our twodimensional model does not represent the flowbandorthogonal (y) components of stress and deformation rates directly, so we parameterise these following Pimentel et al. (2010) assuming no slip at the valley wall: for flow-law coefficient A, longitudinal velocity u, flow-law exponent n, and valley-half width W .By assuming a rectangular glacier cross-section, W is a constant.The strain heating component calculated using this approximation is added to Q str .In the subsequent experiments, the approximate xy strain heating term is small compared to the xz term.
The surface heating term Q m is calculated by assuming that meltwater generation is related only to the difference between the surface air temperature T s and the ice melting temperature T m by means of a constant degree-day factor f dd .The rate of heat capture, per unit height, is where ρ w and L f are the density and latent heat of fusion for water, respectively, and r is a run-off fraction.
The degree-day factor f dd provides a convenient method by which to estimate the summer mass balance based on surface air temperature.The value of the degree-day factor depends to a large extent on the way in which incoming energy is partitioned between different energy balance components (Hock, 2003).Hock (2003) compiled degree-day factors derived for snow at glacierised sites ranging from 2.7 to 11.6 mm d −1 K −1 .Values for ice are typically larger, but are not used here; in our model the degree-day factor is only used to calculate meltwater entrapment (Eq.9) in the accumulation zone, where snow cover is assumed to be perennial.The range chosen (Table 2) spans the commonly reported values tabulated by Hock (2003).
The run-off fraction r allows the removal of a portion of the annual surface melt.The firn captures the remaining meltwater and stores it in a near surface aquifer (cf.Reeh, 1991, for a similar method).Run-off fractions provide a convenient means of estimating internal accumulation, although comparisons with more developed methods find that this approach has limited skill in predicting the thickness of superimposed ice (Wright et al., 2007;Reijmer et al., 2012).Our purpose is to test the influence of a surface heating term that allows for meltwater capture rather than to model melt quantities accurately.The value r = 0.4 has been reported for Greenland near the run-off limit (Braithwaite et al., 1994), but this should vary depending on firn thickness, firn temperature, and summer mass balance.Rabus and Echelmeyer (1998) give estimates of internal accumulation on McCall Glacier that imply high inter-annual variability, although this may be exaggerated by the mercurial accumulation zone conditions on McCall Glacier.The run-off fraction is, therefore, poorly constrained, so with a reference value of r = 0.5, we alter the run-off fraction between 0.2-0.8 in order to evaluate a range of contributions to water entrapment.
The near-surface aquifer is restricted to the accumulation zone and its thickness h aq does not vary annually in the model.This parameter physically represents the thickness of the permeable surface layer through which water can percolate.Due to refreezing and the formation of ice lenses, the near-surface aquifer thickness may be less than the total firn thickness.A suitable choice for the near-surface aquifer thickness depends on climatology.(1994) report a percolation depth of 2-4 m on Greenland, while Fountain (1989) estimates the aquifer thickness on South Cascade Glacier in Washington to be 1.25 m.Firn water in Storglaciären resides in a layer up to 5 m thick, while on Aletschgletscher in Switzerland the firn aquifer is 7 m thick (Jansson et al., 2003).It is reasonable to expect that the aquifer thickness varies spatially, perhaps being thicker at high elevations resulting in a tapered shape.Alternatively, colder temperatures at higher elevations may cause faster refreezing and decrease the thickness of the permeable layer.In light of uncertainties in how to best represent variable nearsurface aquifer thickness, we make the minimal assumption that the near-surface aquifer thickness is invariant in space.
We choose h aq = 3 m as a reference value, and test over the range 0.5-6.0 m.If this assumption is violated, areas where the aquifer is thicker would tend to preserve more liquid water through the winter, while areas where it is thinner would preserve less.This might either reinforce or oppose the gradient in water entrapment implied by melt volumes that decrease with altitude.
In the model, water captured in the near-surface aquifer is stored until it either freezes or exceeds a prescribed drainage threshold ω aq .Above this threshold all water is assumed to contribute to runoff and is removed.The corresponding drainage threshold in the englacial aquifer ω eng (Greve, 1997a) is set to a much lower value to account for a lower porosity in ice compared to the surface layer.In the ablation zone, Q m = 0, as all meltwater is assumed to be removed by the end of the melt season and, therefore, unavailable for entrapment or refreezing.
Englacial water content is poorly constrained by observation, and recent results range from under 1 % (Pettersson et al., 2004) up to several percent (Macheret and Glazovsky, 2000).Although our reference model enforces immediate drainage for water content (ω eng ) above of 1 %, it is likely that the permeability and drainage properties of ice vary spatially, for example as suggested by Lliboutry (1976).Firn porosity has been reported by Fountain (1989) for South Cascade Glacier to be 0.15 with 61 % saturation.We assume that the properties of the near-surface aquifer are similar to that of the firn and use a near-surface maximum water content ω aq = 10 % as a reference value.To quantify model sensitivity to ω aq , we test over a range of 1-15 %.Plausible physical reasons for this variation include variations in accumulation rates and temperature-dependent densification rates.
Within the accumulation zone, surface density is assumed constant in time, and varies with depth according to where ρ f is firn density and C is a constant (cf.Schytt, 1958).
The ablation zone experiences seasonal snow cover, which is represented in the model by changing the density of the upper model layer.We have found the enthalpy difference between an ice column treated in this way and a column with the accumulation and ablation of snow explicitly accounted for to Maximum mass balance 1.5 0.5-2.5 m yr −1 (ice-equivalent) be acceptably small (<0.35K equivalent) for the present purposes.Table 1 lists the model parameters held constant in all simulations.

Ice dynamics
The rheology of ice is described by Glen's flow law which relates strain rate ˙ ij to the deviatoric stress τ ij tensor.The effective stress τ E is the second deviatoric stress invariant.In cold ice, the temperature-dependent flow-law coefficient is computed following the recommendation of Cuffey and Paterson (2010, Ch. 3) as where Q c is the creep activation energy in cold ice, R is the ideal gas constant, T m represents the pressure correction for the melting temperature and T (H ) is the temperature of cold ice as a function of enthalpy.T 0 is an arbitrary reference temperature below which enthalpy takes on negative values.Additional softening associated with non-zero water content may be important in temperate ice.Where there is temperate ice, we multiply the flow-law coefficient by an enhancement factor e w that depends on water content ω (cf.Greve, 1997a), We choose the slope of e w (ω) based on results by Duval (1977) that indicate that A e is roughly tripled with 1 % water content.The tripling of A approximately spans the range of observed strain-rate enhancement in temperate glaciers (Cuffey and Paterson, 2010), so we do not extrapolate further.The applicability of the water-dependent strain-rate enhancement measured in laboratory experiments to glaciers remains an open question that we do not consider here.
Velocities are obtained through vertical integration of the shear strain rate component ˙ xz .We have compared the stress and velocity fields obtained using a "Blatter-type" first-order approximation (FOA) of the momentum balance described by Pimentel et al. (2010) and the zeroth-order shallow ice approximation (SIA).In brief, the SIA reduces the momentum balance in the x-direction to ∂σ xz ∂z = −ρg for vertical shear stress σ xz and ice surface elevation z s .The SIA is most applicable to ice masses with low aspect ratios and small bedrock gradients (e.g.Le Meur et al., 2004).The implications of using the SIA versus the FOA for the present purposes are discussed below.

Boundary conditions
The basal boundary for Eq. ( 3) is treated as a Neumann-type boundary with an enthalpy gradient consistent with The upper boundary condition for Eq. ( 3) is Dirichlet-type where the ice enthalpy is pinned to match either the air temperature or the ice melting point, whichever is lower.We represent annual air temperatures as a sinusoid.A scalar offset T accounts for changes in temperature between model runs.As a function of Julian day t, air temperature T is parameterised as for mean annual temperature T ma at a reference elevation, and vertical lapse rate ∂T /∂z.We ignore shorter period temperature fluctuations because they have a shallower depth of penetration into the ice than the annual cycle.
For the ice dynamics, the surface boundary is a zero-stress boundary, while u = u b at the basal boundary and u b = 0 in most experiments (cf.Le Meur et al., 2004).By imposing u x = 0 at the bed, heat near the bed is advected more slowly and rates of englacial strain heating are slightly elevated compared to experiments in which basal sliding is permitted.A more in-depth investigation of sliding and thermal structure is part of a separate study (Wilson et al., 2013(Wilson et al., , 2012)).

Implementation
Equation (3) steps forward in time on a two-dimensional structured grid that is irregularly-spaced on the vertical axis (z) and regularly-spaced on the horizontal axis (x).The grid spacing in z is finely resolved at the surface and basal boundaries and coarser in the glacier interior.We solve the first term on the right-hand side of Eq. ( 3) using an energy-conserving Crank-Nicolson finite-difference scheme.Because of the small thickness-to-length ratio of glaciers, we omit horizontal diffusion by substituting into the first term on the righthand side This "shallow enthalpy" approximation is identical to that made by Aschwanden et al. (2012) in the Parallel Ice-Sheet Model (PISM).We solve the second term on the right-hand side of Eq. ( 3) in two-dimensions using a flux-limited linear upwind differencing scheme (LeVeque, 1992, Ch. 16).This method is less diffusive than first-order upwind differencing, yet preserves monotonicity in the neighbourhood of large derivatives such as near the ice surface.The model timestep is chosen in the range 30-60 days based on what is empirically found to permit convergence.The timestep in the upwind differencing scheme is permitted to decrease adaptively in order to maintain numerical stability.
Our coupling scheme synchronously steps forward in time in the flow mechanics model and the thermal model.The flow mechanics model computes velocity and stress fields, while the thermal model solves for the enthalpy field.When glacier geometry is permitted to evolve, the ice surface changes based on the mass continuity equation and the prescribed mass balance.Timesteps for the thermal model are limited by the requirement that seasonal changes in Q m (Eq.9) be resolved rather than by stability criteria.We calculate the flowlaw coefficient at every flow-mechanics timestep based on the enthalpy field.

Experimental design and reference models
In order to address the goals discussed above, this study is organised into three experiments.These experiments investigate (1) the primary controls on glacier thermal structure, (2) the sensitivity of thermal structure to environmental and model parameters, and (3) changes in thermal structure accompanying rising air temperatures.As a control for the above experiments, we create two reference models (described below) representing different thermal regimes.In the following section, we outline the methodology behind the three experiments.
We use a simple glacier geometry to isolate the influence of individual environmental and internal variables on thermal structure.Simple glacier geometries also help preserve generality by avoiding effects introduced by irregularities in the prescribed surface and bed topography that might be unique to individual glaciers.We represent the bed as a low-order polynomial function (Fig. 1).
Net balance is approximated as a linear function of ice surface elevation with a prescribed equilibrium line altitude z ELA and balance gradient ∂ ḃ/∂z.We cap the maximum annual accumulation at ḃmax , giving the annual balance function a piecewise-linear shape: where z max = ḃmax (∂ ḃ/∂z) −1 + z ELA .Accumulation is not addressed directly, but rather is implicitly assumed to compensate for the melt calculated in Eq. ( 9) such that the sum matches the prescribed net balance.We experiment with changing the balance gradient (∂ ḃ/∂z) and the balance threshold ( ḃmax ) to simulate glaciers with higher and lower rates of mass turnover.
The steady-state reference models are based on the SIA and incorporate all of the heat sources discussed above.The first reference model (REFT) contains a large volume of temperate ice and arises from the parameter values given in Table 1.The second (REFC) is a colder version of the first produced by shifting the air temperature (T ma ) down by −1.5 K.The REFT model glacier is polythermal (Fig. 1a), with a distribution of temperate ice that is similar to the type "C" configuration illustrated by Blatter and Hutter (1991).Ice within the accumulation zone of the glacier is temperate, while a surface layer of cold ice develops in the ablation area.The bed at the terminus is cold.Similar thermal structure has been observed in Svalbard (Dowdeswell et al., 1984), in Scandinavia (Holmlund and Eriksson, 1989), in the Alps (Vincent et al., 2012;Gilbert et al., 2012), and on the continental side of the Saint Elias Mountains in Yukon, Canada (unpublished data, Simon Fraser University Glaciology Group).
The thermal structure of this reference model is heavily influenced by meltwater entrapment in the accumulation zone.By comparison, meltwater entrapment plays a more limited role in the REFC model (Fig. 1b).Lower surface temperatures decrease the quantity of meltwater production in the accumulation zone, thus, reducing the amount of heat generated at the surface.The REFC model exhibits a type "D" thermal structure as identified by Blatter and Hutter (1991), with a temperate zone in the lower ice column.As in REFT, REFC is frozen to the bed at the terminus.
The SIA neglects lateral and horizontal axial stresses, which is understood to be problematic in mountainous areas where bedrock slopes are large (Le Meur et al., 2004).In regions with steep slopes, this omission leads to unrealistically high deformation gradients and large strain rates.The simple glacier geometry used in this study (Fig. 1) has small (≈ 0.1) bedrock slopes throughout most of the domain, including over the region where modelled ice thickness becomes greatest.This reduces the discrepancy between the SIA and the more correct FOA.Bedrock gradients are steeper than 0.3 over the first 1 km of the model domain, which causes the SIA to deviate more from the FOA in this area.Modelled ice thickness is small in the first 1 km, so the effect of velocity overestimates on strain heating is small.
When the FOA is used to reproduce the REFT and REFC models above (Fig. 1c-d), a similar set of polythermal structures is obtained.Heat content is lower over most of the domain for the FOA model because axial stresses reduce the magnitude of the deviatoric stress tensor.For the REFT model, the FOA ice thickness (Fig. 1c) is greater and the ice extent is shorter by 800 m, consistent with the results of Le Meur et al. (2004).The FOA model exhibits a cold layer that is up to 10 m thicker.The temperate ice fraction in the SIA model is 71 %, while in the FOA model, it is 76 %.The higher temperate ice content in the FOA model owes to the shorter length of the glacier below the equilibrium line altitude and the correspondingly smaller cold ice layer.For the REFC model, the thermal structure is again similar between the SIA and FOA models.The temperate ice fraction in the SIA model is 11 %, while for the FOA model it is 12 %.The length, ice thickness, and heat distribution predicted by the SIA and FOA models, while not identical, are similar enough for the chosen bed geometry that we rely on the SIA models, henceforth.
Although we choose to neglect sliding in this study for the reasons given in Sect.2.1.4,it is desirable to qualify what this omission might imply for our results.To do this, we have added a Weertman-style sliding law (e.g.Cuffey and Paterson, 2010) to the REFT and REFC models, with a sliding coefficient chosen to admit basal velocities in the vicinity of ∼10 m yr −1 .When the ice surface is held fixed, we find that the increase in advection rates caused by sliding thins the cold ablation zone layer in REFT, and causes the basal temperate layer in REFC to be thicker near the terminus.When the surface is permitted to evolve in response to the changing viscosity and material advection rates, the glacier elongates with sliding permitted.The REFC results are similar to above, with the temperate zone thickening near the terminus.In the REFT model, the fraction of temperate ice is lower because rapid sliding rates in the steep upper glacier advects more cold ice.This is not very significant because the actual change in enthalpy is smaller than the large change in temperate fraction would suggest, and the upper glacier remains near the melting temperature.At the broad scales relevant for this study, modest sliding rates would not alter the conclusions presented.

Experiment 1: heat source contributions
In order to investigate the relative importance of different heat sources (Objective 1), we begin with the REFT model and individually remove the contributions from strain heating Q str , meltwater entrapment Q m , and basal heating Q b before recomputing steady-state thermal structure.We test the effect of allowing glacier geometry to evolve in response to changes in ice viscosity governed by Eqs. ( 12) and ( 13), and compare the results to simulations with a fixed glacier geometry.The appropriateness of holding the surface geometry fixed depends on the degree to which thermal structure alters ice fluidity in Eqs. ( 12) and ( 13).Because of the large differences in temperate ice volume in the REFT and REFC models, we examine the effect that flow-coefficient parameterisation has on both.

Experiment 2: parameter sensitivity
Starting from the REFT model, we vary selected parameters in order to explore the sensitivity of steady-state thermal structure to environmental conditions (Objective 2).The parameter ranges given in Table 2 span the spectrum of interesting and physically-meaningful model behaviour.Each parameter is first adjusted independently.To maintain simplicity, we do not consider seasonal perturbations, which may nevertheless be relevant to parameters such as air temperature.We perform tests using both the temperature-dependent flow-law coefficient A (Eq. 12) and the enhanced flow-law coefficient A e (Eq.13).In reality, the parameters in Table 2 are not independent, but considering them as such yields information about the environmental variables controlling thermal structure without complicating the results with multiple causes.Furthermore, considering independent parameters here avoids assumptions about how parameters may be coupled.
Heat flow within glaciers has been described as advectiondominated (characterised by high Péclet numbers) (Aschwanden and Blatter, 2009), but due to the wide range in worldwide glacier velocities, the relative importance of heat transfer by advection compared to diffusion varies.We explore the role played by advection in governing thermal structure by adjusting an additional parameter that is a multiplicative factor (C u ) on advection rate u in Eq. ( 3).We adjust the advection rate directly rather than by changing the mass balance function (Eq.18) or the flow-law coefficient (Eq.12) in order to more clearly isolate experimental variables.Sensitivity tests on C u illustrate the extent to which heat flow in the reference glacier is advection-dominated.They can also be used to investigate the implications of changing flow velocities that result from dynamic behaviour not explicitly considered in the fixed-geometry experiments.
We recognise that some of the variables from Table 2 are correlated.Therefore, we vary air temperature (T ), equilibrium line altitude (z ELA ), and near-surface aquifer thickness (h aq ) together in order to explore the effects of more realistic forcing regimes.To simplify the interpretation of the results, ice geometry is held fixed.
We use the results of the sensitivity tests to draw preliminary conclusions about how glacier thermal structure may evolve in a changing environment.To estimate how nearsurface aquifer thickness and equilibrium line altitude might co-evolve with air temperature, we make use of balance sensitivities: where ḃn is the net balance and ∂ ḃn /∂T can be estimated from field data (e.g. de Woul and Hock, 2005;Oerlemans et al., 2005).Equation ( 19) defines how equilibrium line altitude varies with changes in air temperature.
For Experiments 2 and 3, we introduce the assumption that www.the-cryosphere.net/7/167/2013/The Cryosphere, 7, 167-182, 2013 the near-surface aquifer thickness is related to net balance closely enough that net balance sensitivity estimates are applicable and equivalent to aquifer thickness sensitivities: We recognise that this simple parameterisation is an important assumption, but it qualitatively captures the behaviour that we hypothesise.Namely, we expect that as temperatures rise the near-surface aquifer will grow thinner and this rate of change will to a first approximation be proportional to the change in accumulation rates.The percolation pathways within the aquifer also involve ice lenses formed by seasonal refreezing of meltwater (Jansson et al., 2003).We do not consider lateral transport within the near-surface aquifer, but we speculate that the presence of ice lenses will have a limiting effect on aquifer thickness that brings the aquifer sensitivity into closer agreement with the annual balance sensitivity (e.g.Eq. 20).
Where glacier geometry is held fixed in Experiment 2, the results of the sensitivity tests do not represent physically consistent thermal regimes and are not directly representative of future thermal structures.All sensitivity tests are, therefore, repeated with a freely evolving ice surface.

Experiment 3: prognostic modelling
Transient feedbacks between variables such as mean annual air temperature and accumulation zone extent can be expected, but are not represented in the experiments above.In order to capture such feedbacks and make realistic projections of thermal structure (Objective 3), glacier geometry must be permitted to evolve.We perform prognostic simulations with a range of transient climate forcing scenarios.These scenarios are distinguished by the extent to which the winter balance offsets increasing summer ablation.The initial conditions are the REFT and REFC models.We prescribe an average annual air temperature that increases linearly by 2.5 K over 100 yr and then stabilises.For each model timestep, the near-surface aquifer thickness and equilibrium line altitude are adjusted to track the prescribed temperature according to Eqs. ( 19) and (20).A fraction of the ablation response to changing temperature is assumed to be offset by changing winter balance.We text 19 possible winter balance responses to warming that offset ablation by fractions spanning 5 %-95 %.

Evaluation metric
To make quantitative comparisons between different simulations, we use two metrics to describe the modelled glacier thermal structure: (1) equivalent temperature difference relative to a given reference model and (2) temperate ice fraction.The former converts the enthalpy difference between two models into an equivalent temperature field (in Kelvin): This is useful for comparing experiments in which the glacier geometries are identical.The second metric is a simple area fraction of temperate ice along the modelled flowband.Where applicable, both metrics are used.

Experiment 1: heat source contributions
With the temperature-dependent flow-law coefficient A (Eq. 12), the modelled enthalpy without strain heating (Q str = 0) is smaller relative to the reference run (Fig. 2c).With the enhanced flow-law coefficient A e (Eq.13), this difference is slightly larger.The temperate fraction and mean equivalent temperature difference over the entire flowband domain are presented in Table 3.
Even in the upper half of the ice column where stresses and deformation rates are low, equivalent temperatures are lower in general when strain heating is neglected (Fig. 2c).Lower deformation rates at depth lead to lower flow velocities in the upper ice column.Lower velocities reduce heat advection from the accumulation zone source Q m derived from meltwater entrapment.This effect exists whether using flowlaw coefficient A or A e , but is greater with the latter because viscosity becomes a function of water content as well as temperature.
The omission of meltwater entrapment (Q m = 0) causes a large change in the modelled thermal structure (Fig. 2b).The resulting distribution of temperate ice is similar to the REFC model, in which meltwater entrapment has been physically reduced by lowering surface temperatures.The large mass of temperate ice in Fig. 1a becomes limited to the deepest parts of the glacier ablation zone, and the bulk of the ice remains cold.The temperate ice distribution is most similar to the type "D" structure in Blatter and Hutter (1991).In the accumulation zone, the near-surface equivalent temperature is much lower than in the reference run (Fig. 2d).The equivalent temperature differences are again slightly larger using the enhanced flow coefficient A e .Because of the small amount of temperate ice, the enhancement factor in Eq. ( 13) plays a small role.
Relative to strain heating and meltwater entrapment, basal heat sources (not shown) play only a small role.Sufficient basal heating causes the bed to reach the melting point.Temperate conditions do not extend to the glacier interior because of the negligible thermal diffusivity κ t .The insensitivity of the thermal structure to basal heat sources partially justifies our choice not to explicitly include frictional heating from basal sliding.Our simplified model does not capture changes in sliding that might occur in concert with changing basal temperatures.
Allowing the glacier geometry to evolve does not have a strong influence on the results of this experiment.In the absence of strain heating, the cold-temperate transition in the REFT model becomes deeper compared to the control near the glacier terminus.In the absence of meltwater entrapment, the glacier grows thicker and longer due to the higher viscosity of cold ice.A thin temperate zone forms at depth in the lower half of the glacier.This temperate layer is thicker than in the fixed geometry model because of the higher stresses in the larger steady-state glacier (Fig. 3a).
The effect of water content on ice rheology as represented by A e has a large impact on glaciers with a high temperate ice fraction, and a plausibly small impact on those that are mostly cold.Unlike other sources of strain enhancement such as lattice preferred orientation and ice impurity content, water content is directly connected to thermal structure, and by extension, to climate conditions.When glacier geometry evolves freely, both REFT and REFC models with the enhanced flow-law coefficient A e exhibit a steady-state that is thinner than that using A because of the higher fluidity in the large regions of temperate ice.Using the enhanced flow-law coefficient A e with the REFC model yields a glacier thickness 3 % (8 m) smaller than with A (not shown).The thickness of the REFC basal temperate layer, on the other hand, decreases by 40 % (14 m), such that simulations with A predict a thicker temperate layer (Fig. 3a).This thinning with the enhanced coupling of A e is consistent with the findings of Aschwanden et al. (2012).The cold layer thickness in the REFT model, which depends on both advection rates and strain heating, is slightly greater with A e than with A (Fig. 3b).

Independent variables
Changes in boundary heat sources/sinks and internal heat generation require the overall thermal regime to shift in response.Thermal structure within the control models also varies slightly with the choice of A or A e .The sensitivities of the steady-state control models to changing environmental parameters is similar regardless of whether A or A e is chosen, so in the following section we focus on results based on A alone (Figs. 4,5 and 6).
Shifts in air temperature exert a strong control on thermal structure (Fig. 4a).The transition from a fully cold to a mostly temperate glacier occurs over an air temperature range of approximately 3 K.An intermediate thermal structure develops between the two end-member models with a temperate core heated partly by strain (Fig. 7a and b).With further warming, meltwater entrapment in the lower accumulation zone produces temperate ice through the full thickness of the glacier in the central region of the flowband (Fig. 7ce).Cold ice advected from high elevations persists as a cold region upglacier from the temperate zone.
The expansion of the temperate ice region with increasing air temperature occurs by two mechanisms.First, in a warmer climate less heat is lost during the winter and the cold layer that forms in the ablation zone does not penetrate as deeply into the ice.Secondly, larger amounts of heat derived from meltwater entrapment are added to the glacier in the accumulation zone and advected downstream.Within the model, this second effect is partially muted when the firn aquifer becomes saturated, however, the shortening of the cold season associated with higher T causes incrementally less heat to be lost to the atmosphere over the entire elevation range.
In Fig. 4a, there is a sharp transition from high to low temperate ice fractions at air temperatures below that used to produce the reference model (REFT).In REFT, upstream heating from meltwater entrapment is the source of much of the temperate ice in the glacier interior, so eliminating this heat source produces a transition to a glacier with a small temperate ice fraction.The transition is not as rapid if defined in terms of the metric mean K , indicating that englacial heat storage is not strongly affected by the transition from less diffusive temperate ice to more diffusive cold ice.Although the Péclet number typically drops as ice cools to submelting temperatures, the effect on total heat storage within the glacier is small.
The thickness of the near-surface aquifer h aq (Fig. 4b), the degree day factor f dd (Fig. 4c), and the run-off fraction r (Fig. 4e) are important for similar reasons as surface temperature.In the case of a thin surface aquifer, less water from the previous melt season is entrapped, and refreezing and heat loss to the atmosphere occur more efficiently during the cold season (Fig. 7f-i).A thicker aquifer captures more water and preserves more energy at depth because of the insulating properties of the overlying snow and firn (Fig. 7k).The thermal structure of REFT is insensitive to near-surface aquifer water content thresholds ω aq ≥ 5 % (Fig. 4g) and declines below that.The degree-day factor and run-off fraction are not independent parameters as implemented in Eq. ( 9), and together control the amount of meltwater available for entrapment.Reducing the extent of the accumulation zone by raising the equilibrium line (Fig. 4d) has an effect on temperate ice generation because it diminishes the region over which heat can be added.
There is a tendency for temperate ice generation in the accumulation zone to be greatest at intermediate elevations.In this region, temperatures are high enough in the summer to produce significant quantities of melt, and burial rates are high enough to advect large amounts of heat into the glacier before it is lost to cooling in the winter.Nearer the equilibrium line, submergence rates are lower and the ratio of vertical advection to diffusion is smaller.Our assumption of a constant near-surface aquifer thickness causes the upper glacier transition from complete refreezing to producing temperate ice to be different than it would be in the case of a tapered aquifer.Therefore, different aquifer geometries could conceivably alter the zone in which heat from meltwater entrapment is pumped into the glacier.
Increasing the maximum permitted water content ω eng results in higher fractions of temperate ice (Fig. 4f).The ablation zone cold layer thins due to the higher volume of water available for refreezing at the cold-temperate transition surface.The associated higher heat flux requires a steeper thermal gradient such that the cold-temperate transition is nearer the glacier surface.The increase in temperate ice fraction begins to level off with a water content threshold ω eng of 2-3 %.The mean equivalent temperature difference relative to the REFT model increases steadily with ω eng because more heat is stored within the glacier as liquid water.In contrast, basal heat flux (Q b , not shown) does not have a significant effect.This is consistent with the results from Experiment 1.
Additionally, we explore the effect of altering the rate of heat advection by a constant coefficient C u across the glacier.In assigning C u = 1, the velocity field used for energy advection is no longer physically consistent with the glacier geometry.Nevertheless, this experiment is useful for illustrating the effect of varying velocity regimes on thermal structure.a surge), advection will cause the thermal regime to move toward the results developed in these simulations.For low advection rates, the ablation zone cold layer penetrates deeper into the glacier, restricting the extent of temperate ice derived from the accumulation zone.With high advection rates, the resulting thermal structure is similar to that with high allowable water content (Fig. 4h).In both cases, the rate of water transport to the cold-temperate transition increases, causing the transition to occur nearer the glacier surface.
In a final pair of single-parameter sensitivity tests, we investigate the effect of adjusting the vertical mass balance gradient (∂ ḃ/∂z) and the maximum balance threshold ( ḃmax ).In these experiments, the differences in glacier geometry are sometimes large enough that we report results with a freely-evolving ice surface (Fig. 5).When the balance gradient (∂ ḃ/∂z) is small, mass turnover within the glacier is low.The corresponding lower advection rates cause temperate ice to be largely constrained to the upper glacier, and the area of the modelled flowband that is cold is large (Fig. 5a).As the balance gradient rises, the cold ice area drops slightly and the temperate ice volume rises steeply.The mass balance threshold ( ḃmax ) affects the thermal structure largely by restricting glacier accumulation.At low values, the glacier is thinner and flows more slowly, which causes the temperate area of the flowband to be small relative to the REFT control model (Fig. 5b).As the balance threshold rises, the cold ice area stays nearly constant, but the temperate ice area rises.
From the combined results of Experiments 1 and 2, we find that in cold climates and environments where meltwater is not efficiently captured in the accumulation area, strain heating represents a primary control on temperate englacial zones.This heat source is greatest in the deeper part of the ice column.If surface ablation rates are high enough, the layer of ice warmed by strain-heating will eventually be near enough the surface to lose much of this heat to the atmosphere.Such a layered structure (Fig. 7a-b and f-g) is similar to thermal structures observed in Svalbard (Björnsson et al., 1996), the Canadian Arctic (Blatter and Kappenberger, 1988), and in ice streams draining large continental ice sheets (Truffer and Echelmeyer, 2003).Alternatively, when environmental conditions permit meltwater entrapment at the surface, latent heat quickly becomes a dominant heat source.This situation (Fig. 7d-e and i-j) is similar to that observed in glaciers with Each line in (a) represents a scenario with a different net accumulation response to a single function for temperature (b).Near-surface aquifer thickness h aq decreases and equilibrium line altitude z ELA rises when net balance is lowered (all scenarios).Models in which winter accumulation offsets 20 %, 50 % and 80 % of the increased summer melt are shown by the solid, dashed and dotted bold lines, respectively.For reasons described in the text, the lines are terminated when glacier length falls below 3 km.large temperate ice zones, such as Storglaciären (Pettersson et al., 2004).
Model sensitivity to changes in variables depends on the choice of reference glacier (e.g.REFT versus REFC).In terms of temperate fraction, the high sensitivity of meltwaterdominated polythermal glaciers (such as REFT) to small perturbations develops from the large amounts of heat potentially captured (Q m ) or lost through the glacier surface.

Coupled variables
The previous results demonstrate that in the absence of other changes, higher air temperatures may cause polythermal glaciers to become more temperate.Reductions in firn thickness and accumulation area extent have the opposite effect (Fig. 6).In the case of a thick near surface aquifer (Fig. 6a) and a large accumulation area (Fig. 6b), surface air temperature acts as a nearly independent control on thermal structure primarily through the effect of meltwater entrapment.The situation reverses when aquifer thickness becomes less than ∼2 m or when the accumulation area ratio is less than about one-third.
Slices from the parameter space between surface air temperature, near-surface aquifer thickness and accumulation area have been mapped in Fig. 6.These slices show how steady state glaciers in various parts of the parameter space will tend to evolve as conditions change.There are many plausible trajectories across the parameter space, however, a set based on Eqs. ( 19) and ( 20) has been mapped as dashed lines.In drawing these lines in Fig. 6a, we make the assumptions that near-surface aquifer thickness is related to the annual net balance (in firn-equivalent) and that changes in net balance can be represented by a scalar mass balance sensitivity (Eq.19).Published mass balance sensitivities over a 1 K range vary widely, so we choose ∂b n /∂T = 0.5 m K −1 yr −1 (w.e.).For a warming environment, the quasi-steady-state trajectories through the h aq -T and AAR-T parameter spaces are non-monotonic (Fig. 6a and b), in contrast to those through the AAR-h aq parameter space.If a glacier is in a region of the T parameter space where it is insensitive to changes in air temperature (i.e.T < −2.0 for REFT) only the correlated changes in accumulation area and near-surface aquifer thickness are important.In this case, increasing air temperature will produce a reduction in the fraction of temperate ice within the glacier (similar to Fig. 6c as h aq and AAR are reduced).Alternatively, in a scenario where the accumulation area ratio is roughly static, but the firn thinning and temperature rise occur (as in Fig. 6a), a glacier may become more temperate before cooling again as meltwater entrapment is further inhibited.
The sensitivity tests provide a preliminary estimate of how thermal structure may respond to changing climates (Objective 3).A 1 K increase in temperature for the REFT model corresponds to a roughly 8 % increase in the amount of temperate ice in the absence of changes to any other parameters (Fig. 4a).At the same time, if rising temperatures increase summer ablation, the near-surface aquifer thickness (Fig. 4b) should decrease and the equilibrium altitude should rise (Fig. 4d).Assuming that net balance sensitivity to temperature ∂b n /∂T = 0.5 m K −1 yr −1 , the combined effects of aquifer thinning and accumulation zone reduction sum to a much larger decrease in temperate ice fraction than the increase due to increased air temperature alone.This calculation may be altered if a higher winter balance accompanies rising temperatures, so the potential exists for polythermal glaciers to become either colder or warmer in a warming environment (cf.Rippin et al., 2011).

Experiment 3: prognostic modelling
The final experiment examines the transient responses of the REFT and REFC models to changing climate.Equilibrium line altitude z ELA co-varies with a prescribed air temperature evolution in a manner consistent with Eq. ( 19).With rising air temperature, the model net balance decreases linearly.Near-surface aquifer thickness h aq either changes according to Eq. ( 20) or is held fixed.Hypothetical increases in winter www.the-cryosphere.net/7/167/2013/The Cryosphere, 7, 167-182, 2013 balance are prescribed to offset ablation predicted by Eq. ( 19) such that a high winter balance diminishes the effect of rising air temperature on z ELA and h aq .
For the REFT model, a wide range of thermal responses is possible in a warming climate (Fig. 8).Many trajectories show decreasing temperate fraction over time, with this effect being most pronounced when winter balance does little to offset increased summer melt.The glacier length also falls.In shorter glaciers, a large portion of the bed has a high slope that the SIA is ill-equipped to handle.Furthermore, the glacier length relative to the fixed horizontal discretisation becomes small.Therefore, we do not include results from glaciers <3 km long in our analysis.When large increases in winter balance are prescribed, the fraction of temperate ice (as well as the glacier volume) remains relatively steady.More than 80 % of the ablation increase must be offset by increased accumulation in order to maintain or increase the temperate ice fraction for REFT.In this scenario, an 80 % accumulation offset for 2.5 K warming is equivalent to an increase in winter balance of 1 m [w.e.].
The thermal evolution of the colder REFC model is different, with many of the high accumulation offset scenarios resulting in increased temperate ice fraction.Rising air temperatures increase the meltwater availability in the accumulation area.The increased meltwater is captured in scenarios where the near-surface aquifer thickness and extent remain large.For smaller winter accumulation sensitivities (balance offsets less than ∼50 %), the glacier remains dominated by cold ice.In many of the scenarios where winter balance increases little, both REFT and REFC models become small and thin.
In order to test the impact of the assumption in Eq. ( 20), we perform the experiment above with the alternative assumption that the aquifer thickness (h aq ) is constant.The resulting thermal evolution for REFT is similar to above, however, the temperate ice loss is milder.In this case, the effect of the rising equilibrium line (z ELA ) is sufficient to cause a significant reduction in temperate ice extent.For the REFC model, the results are slightly different; a number of balance scenarios result in rising temperate ice fractions as more meltwater is captured in the near-surface aquifer.Therefore, there is uncertainty because the future evolution of thermal structure depends on the behaviour of the near-surface aquifer.Although ∂h aq /∂T is poorly-known, we consider the relationship in Eq. ( 20) to be more realistic than assuming no change.
Based on an assumed accumulation increase of 10 % and a temperature rise of 1 K, de Woul and Hock (2005) find that increased winter balance offsets increased ablation by approximately 20 %.The range of accumulation offsets for glaciers north of 60 • N that they report stretches from 54.4 % to less than 5 %.Our results suggest that offsets > 85 % for the REFT model and > 60 % for the REFC model are required to maintain or increase the temperate ice fractions in the prescribed warming scenario (2.5 K warming).This leads us to conclude that many polythermal alpine glaciers will experience a net cooling as climate warms.

Conclusions
Simple two-dimensional numerical experiments with synthetic glacier geometries reproduce a variety of thermal regimes that have been observed in existing glaciers.These structural configurations exist along a realistic range of environmental parameters.
For a small polythermal reference model dominated by temperate ice (REFT), meltwater entrapment comprises the primary source of heat to the glacier.Meltwater entrapment provides a mechanism through which changes in surface conditions have a direct influence on thermal structure.Strain heating plays a less important role for the glacier geometries tested here.Basal heat fluxes are a relatively unimportant control on temperate ice volume, but they have an effect on basal temperatures.
In the polythermal regime embodied by REFT, environmental factors affecting the formation and entrapment of meltwater in the accumulation zone strongly affect thermal structure.The quantity of surface melting, the equilibrium line altitude and the near-surface aquifer thickness have competing roles in controlling the quantity of temperate ice produced.Thinning firn layers and retreating accumulation zone extent have the potential to cut off the supply of latent heat to the glacier, thereby prompting transitions in thermal regime to reduced or nonexistent temperate zones.The importance of meltwater entrapment in driving thermal structure makes better constraints on accumulation zone parameters desirable.
If increases in accumulation partially offset increases in ablation with rising temperature, the loss of temperate ice volume slows.The increase in accumulation required for a REFT-or REFC-type glacier to preserve or enlarge temperate ice zones is likely beyond that predicted for many glacierised environments.For this reason, many polythermal mountain glaciers are likely to experience cooling in warming climates, particularly where the thermal regime is presently dominated by meltwater entrapment.
In polythermal glaciers with a small volume of temperate ice, the rheological softening of ice due to non-zero water content has only a limited effect on modelled geometry and thermal structure.Where large temperate regions and high water contents exist, rheological softening leads to substantial differences in simulated geometry and thermal structure.These predicted differences for temperate ice masses highlight uncertainties surrounding the formulation of the flowlaw coefficient in real glaciers.

Fig. 2 .Fig. 3 .Fig. 2 .
Fig. 2. Distributions of enthalpy and equivalent temperature difference (∆K ) for REFT with and without strain heating Qstr (a,c) and meltwater entrapment Qm (b,d).Model runs use flow-law coefficient A as in Eq. (12).The dashed lines in (a,b) denote the cold-temperate transition.Darker areas in (c,d) indicate the largest differences in ice enthalpy relative to the REFT reference model.Enthalpy difference in (c,d) in equivalent temperature according to (Eq.21).

Fig. 2 .Fig. 3 .
Fig. 2. Distributions of enthalpy and equivalent temperature difference (∆K ) for REFT with and strain heating Q str (a,c) and meltwater entrapment Q m (b,d).Model runs use flow-law coefficient A a (12).The dashed lines in (a,b) denote the cold-temperate transition.Darker areas in (c,d) indicate th differences in ice enthalpy relative to the REFT reference model.Enthalpy difference in (c,d) in eq temperature according to (Eq.21).

Fig. 3 .Fig. 4 . 28 Fig. 4 .
Fig. 3. Effect of flow-law coefficient parameterisation on glacier thermal structure.Glacier geometry evolves in these simulations.The basal temperate layer thickness (a) is normalised to ice thickness.The cold layer thickness in (b) is equivalent to the coldtemperate transition depth.Note that the x-axes begin near the middle of the glacier in order to focus on the cold ablation zone layer.

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 5 .
Fig. 5.The results of varying (a) mass balance gradient and (b) maximum mass balance, given in terms of mean equivalent temperature (solid line, dots) temperate fraction (dashed line, crosses), and flowband area (grey line).The glacier geometry is permitted to change in response to the changing net balance profile.The values used for the reference model are indicated by the vertical grey bars.

Fig. 6 .
Fig. 6.The effect of varying pairs of variables on temperate ice fraction (contour interval is 0.1).Equilibrium line altitude has been recast as accumulation area ratio.Air temperature offset ( T ) and aquifer thickness (h aq ) are co-varied in (a), air temperature offset ( T ) and AAR in (b) and aquifer thickness (h aq ) and AAR in (c).The simulations shown use the enhanced flow coefficient (A e ).The circles indicate the reference parameter combinations.Dashed lines denote hypothetical trajectories through the parameter space based on a linear mass balance sensitivity and a constant lapse rate.See text for details.

Fig. 6 .Fig. 7 .Fig. 7 .
Fig. 6.The effect of varying pairs of variables on temperate ice fraction (contour interval is 0.1).Equilibrium line altitude has been recast as accumulation area ratio.Air temperature offset (∆T ) and aquifer thickness (haq) are co-varied in (a), air temperature offset (∆T ) and AAR in (b), and aquifer thickness (haq) and AAR in (c).The simulations shown use the enhanced flow coefficient (Ae).The circles indicate the reference parameter combinations.Dashed lines denote hypothetical trajectories through the parameter space based on a linear mass balance sensitivity and a constant lapse rate.See text for details.

Fig. 8 .
Fig. 8. Temporal evolution of the fraction of temperate ice with changing environmental conditions for time-dependent models.Each line in (a) represents a scenario with a different net accumulation response to a single function for temperature (b).Near-surface aquifer thickness h aq decreases and equilibrium line altitude z ELA rises when net balance is lowered (all scenarios).Models in which winter accumulation offsets 20 %, 50 % and 80 % of the increased summer melt are shown by the solid, dashed and dotted bold lines, respectively.For reasons described in the text, the lines are terminated when glacier length falls below 3 km.

Table 1 .
Physical constants and model parameters.

Table 2 .
Environmental parameters varied in model sensitivity tests.

Table 3 .
Results of heat source removal (Experiment 1) with fixed glacier geometry and flow-law coefficient A (Eq. 12).