Skip to main content Accessibility help


  • Access
  • Cited by 13


      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Sensitivity of the frozen/melted basal boundary to perturbations of basal traction and geothermal heat flux: Isunnguata Sermia, western Greenland
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Sensitivity of the frozen/melted basal boundary to perturbations of basal traction and geothermal heat flux: Isunnguata Sermia, western Greenland
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Sensitivity of the frozen/melted basal boundary to perturbations of basal traction and geothermal heat flux: Isunnguata Sermia, western Greenland
        Available formats
Export citation


A full-stress, thermomechanically coupled, numerical model is used to explore the interaction between basal thermal conditions and motion of a terrestrially terminating section of the west Greenland ice sheet. The model domain is a two-dimensional flowline profile extending from the ice divide to the margin. We use data-assimilation techniques based on the adjoint model in order to optimize the basal traction field, minimizing the difference between modeled and observed surface velocities. We monitor the sensitivity of the frozen/melted boundary (FMB) to changes in prescribed geothermal heat flux and sliding speed by applying perturbations to each of these parameters. The FMB shows sensitivity to the prescribed geothermal heat flux below an upper threshold where a maximum portion of the bed is already melted. The position of the FMB is insensitive to perturbations applied to the basal traction field. This insensitivity is due to the short distances over which longitudinal stresses act in an ice sheet.


At geologic timescales, high-latitude countries in the Northern Hemisphere will likely experience future glaciations. The long-term storage of nuclear waste in deep geologic repositories can potentially be impacted by a glaciation via the ice sheet’s influence on the subglacial and proglacial groundwater system. It is therefore important to consider subglacial hydrological processes and the role ice sheets play in driving groundwater systems, when designing safe storage systems in northern locations. Subglacial hydrological processes become active and recharge the groundwater system only where the bed of an ice sheet is melted. Understanding the spatial pattern of thermal conditions of an ice sheet’s bed is therefore an important design criterion for responsible nuclear waste disposal in countries such as Sweden, Finland and Canada.

The Greenland ice sheet (GrIS) is a present-day analog to future ice sheets in Scandinavia and Canada. The thermal state of the bed of GrIS and the accumulation of subglacial water has been investigated by a variety of methods, but remains poorly constrained. Direct observations via drilling show that melted conditions exist near the western margin (Lüthi and others, 2002), as well as at a north-central location near the ice-sheet divide (Andersen and others, 2004). Conversely, frozen conditions have been noted at point locations spanning the ice sheet from Camp Century near the northwestmargin (Weertman, 1968), to the centrally located Greenland Icecore Project (GRIP) core and the southeast Dye 3 core (Dahl-Jensen and others, 1998). The spatial extent of melted bed conditions, as determined by the few point observations, has been extended via interpretation of ice-penetrating radar. Fahnestock and others (2001) derived spatially variable basal melt rates exceeding 0.15 ma–1 in central GrIS through interpretation of internal radar layering. Using bed reflectivity power as a proxy for basal water content, Oswald and Gogineni (2008) suggested a spatially heterogeneous basal water distribution along radar transects of GrIS.

Spatially comprehensive estimates of basal conditions are offered by ice-sheet model output. Greve and Hutter (1995) investigated the sensitivity of the basal temperature field on GRIS to variations in a uniform geothermal heat flux. Their results suggest that, while increasing heat flux caused an inland migration of temperate basal conditions, the interior remained frozen, even under the highest heat-flux scenario (54.6mWm–2). This was complemented by a follow-up investigation of the basal temperature field, matching geothermal heat flux to point observations, which implied that the majority of the ice-sheet bed was at the pressure-melting point (Greve, 2005). In addition to geothermal heat flux, the sensitivity of GrIS basal conditions to changes in surface temperatures and mass balance was investigated by Huybrechts and others (1996), who found basal conditions show a pronounced sensitivity to steady-state changes in temperature and mass balance: for example, a 10 C drop in surface temperature resulted in a freezing of the majority of the ice-sheet bed; however, the drop in surface mass balance associated with a 10 C lowering of surface temperature resulted in temperate conditions over nearly 60% of the bed. Transient simulations over the last two glacial cycles show most of GrIS exhibits frozen conditions at the bed at some point in time. While the models employed by both Greve and Huybrechts were three-dimensional, both were mechanically limited to the shallow-ice approximation.

In summary, previous studies suggest a spatially distinct frozen/melted boundary (FMB). The location of the FMB at the bed is the result of a balance between heat sources concentrated near the bed (frictional heat from sliding, geothermal heat flux and strain heating) and the introduction of colder ice through diffusive and advective processes. In the present study, we investigate the sensitivity of the FMB not only to geothermal heat flux, but also to changes in cold ice advection resulting from ice motion, including basal sliding. Sensitivity is investigated with a steady-state, thermomechanically coupled, two-dimensional flowline model which solves the full-stress equations (a vertically explicit solution that includes membrane stresses). We apply this model to a profile of Isunnguata Sermia (Fig. 1), a terrestrially terminating glacier in western GrIS. The model is brought into agreement with observation by using adjoint methods for evaluating gradients of an objective function. Motivation for selecting a terrestrially terminating glacier stems from the fact that the majority of GrIS is land-terminating, and such a profile removes additional physical complexities relating to marine-terminating ice. Using the steady-state glacier geometry and surface velocity field, we examine the interactions of heat sources that dictate the stability of the FMB under different assumptions about geothermal heat flux and the basal traction fields.

Fig. 1. Study site displaying model profile (red line) from the ice-sheet divide through Isunnguata Sermia to the western margin. Surface elevation contours (blue lines) are given in meters above sea level, and interpolated from Bamber and others (2001). The yellow contour at 1500ma.s.l. represents the approximate equilibrium-line altitude (ELA), according to Van de Wal and others (2008).


Field equations

Our model is built upon the continuum mechanical formulation of the laws of conservation of mass, momentum and energy for an incompressible fluid. These are, respectively,




where u is the velocity vector, σ the stress tensor, θ the temperature and Φ sources of heat generation in the ice. Physical constants cp , k i, ρ and g are defined in Table 1. Analysis is restricted to the x-z plane, or the vertical profile, making , where î and are unit vectors in the x and z directions, respectively.

Table 1. Parameters and physical constants used in the model

Conservation of momentum and mass

The constitutive relation for ice takes the form


where τij is the ij element of the deviatoric stress tensor, which is defined by τij = σij – pδij , with δij the Kronecker delta function. Isotropic pressure is defined as p = –⅓∑ i σii . є ij represents the corresponding element of the strain-rate tensor and η the viscosity. The strain-rate tensor is given by, and related to, velocity gradients as follows:


A non-Newtonian rheology is used for ice:


with , or the second invariant of the strain-rate tensor, and a regularization parameter introduced to avoid a singularity at zero strain rate. Glen’s flow law (Paterson, 1994) gives n = 3. A(θ* ) is the flow law rate factor, given by Paterson and Budd (1982):


where θ* is the homologous temperature, defined by θ* = θ + bp, and R is the universal gas constant.

Under the assumption of steady state, the velocity of the ice is then determined from Stokes flow confined to the x-z plane,


and the conservation of mass, ∇u = 0.

Conservation of energy

Φ, the term in Equation (3) which represents internal heat generation, is computed as


Under the assumption of steady state and uniform thermal conductivity, the temperature of the ice is given by


Boundary conditions

Boundary conditions are applied to three distinct regions on the boundary of Isunnguata Sermia: (1) the surface; (2) the bed; and (3) a vertical boundary at the ice divide.

Conservation of momentum and mass boundary conditions

The surface of the glacier upholds the neutral or stress-free boundary condition:


where is the outward normal unit vector.

The bed of the glacier is subjected to a Weertman-style sliding law, where basal velocity and shear stress are related as


where β 2 is a positive scalar, spatially variable parameter representing the magnitude of frictional forces at the bed and τ b is given by


evaluated at the base of the glacier. We constrain the sliding velocity to be tangential to the bed, that is u =0.

The vertical boundary at the divide is subject to a symmetry boundary condition:



where is the unit vector tangent to the divide.

Conservation of energy boundary conditions

The bed of the glacier is subject to an inward heat flux given by


where q g is the geothermal heat flux, q f is heating due to sliding friction and q l is latent heat associated with the melting of ice. q g is taken as 42 mWm–2, unless otherwise stated. Frictional heat is calculated as


Latent heat is given by


This heat interacts with the ice via the Neumann boundary condition


Note that the inclusion of the latent heat term serves as a temperature constraint on the ice by counteracting the inward flux from geothermal heat and frictional heat when the basal ice is at the pressure-melting point.

The surface temperature of the glacier is inferred from the dataset of Ettema and others (2009), and is imposed as a Dirichlet boundary condition. The vertical boundary at the divide is thermally insulated such that [−k i ∇θ] = 0.

Model domain

The geometry for the model domain was derived from the surface elevation and thickness data of Bamber and others (2001). Since the model used here considers only a vertical profile, we selected a streamline from the surface velocity data presented by Joughin and others (2010). We employed cubic splines to interpolate the glacier geometry between data points, which were spaced at 5 km.

Due to the discrete nature of the original dataset, the profile surface contained numerous artifacts, manifested as irregularities in slope. In order to produce a more reasonable surface, we implemented a free-surface evolution scheme, and allowed the model geometry to relax for 50 years. The high driving stresses associated with the slope irregularities quickly diffused, yielding a surface free from the original artifacts, but still consistent with the data and model physics.

Numerical considerations

The model uses the finite-element method to solve the field equations subject to the boundary conditions. Lagrange quadratic elements are used (Hughes, 2000), allowing second derivatives of the velocity to be computed accurately. The nonlinearity resulting from the viscosity (Equation (6)) is resolved using the modified Newton’s method iterative solver (Deuflhard, 1974). The resulting linear systems were solved with UMFPACK (Davis, 2004). Model-specific parameters are summarized in Table 2. All numerical work was carried out in the Comsol Multiphysics modeling environment, a commercial package for finite-element analysis of general partial differential equations.

Table 2. Quantities of importance for model numerics

Modeling assumptions

Several assumptions were made in the development of this model, and results must be understood with these in mind. The assumptions are as follows:

The datasets used in the generation of the model domain geometry are sufficiently accurate, and the surface smoothing used to reduce artifacts does not introduce additional errors larger than those resulting from artifacts in surface geometry.

Stresses acting transverse to the dominant flow direction are small. This is necessary due to the effect that these stresses, and associated strains, have on the rheologic properties of the ice. Given the profile’s location at the center of the ice catchment, and the uniform width of the streaming feature, this assumption is likely to be valid.

The steady-state solution generated by the data-assimilation process is a reasonable representation of a long-term configuration for the model domain. This assumes that the modeled region of GrIS was not in a transient state when the data were collected.

A constant geothermal heat flux is an appropriate parameterization of the real phenomenon across the modeled domain. This is to say that, given the spatial scale under consideration, variability in geothermal heat flux is either of a sufficiently low resolution to be considered in an average sense, or of a sufficiently large scale that it is essentially constant.

Steady-state solutions which include the data-assimilation process are sufficient for probing the sensitivities of the system with respect to changes in the basal boundary. A more complete treatment would entail the evolution of the free surface to determine the ultimate outcome of the perturbation, but that is beyond the scope of this work.

Data assimilation and model initialization

When modeling ice dynamics, there are two issues that must be addressed before numerical experiments can be conducted. Firstly, fields which have not been directly measured but are significant in computing flow must be estimated. For instance, the internal distribution of temperatures is critical to ice dynamics, but is at best known at a few boreholes. We refer to this process as ‘model initialization’. Secondly, the initialized model should be in agreement with measurements that are available. We refer to this as ‘data assimilation’.

Our strategy in this paper is to use steady-state solutions to conservation equations to initialize the model, subject to the constraints introduced by the data-assimilation process. This is not a new idea: MacAyeal (1993) introduced control methods in the context of ice-sheet modeling. Here, we extend the concepts to solutions which incorporate the full flowline stress balance.

For data assimilation, we use the adjoint of the linear operator to compute derivatives of an objective function, and use these slopes to minimize the function. We have defined an objective function in terms of difference between the observed (u obs(x)) and modeled (u mod(x)) surface velocities,


which will be differentiated with respect to a parameter, p, that we vary in order to minimize the objective function. In this case the parameter is p = β(x)2, or the basal traction. Our introduction follows that of Strang (2007, p. 678–684).

‘Chain rule’ differentiation yields


where u is a solution vector containing both velocities and temperatures. The key to efficient calculation of the derivatives of the objective function is writing


or, recognizing that the objective function is linear in u. It is now possible to write the gradient as


where c T A 1 is the result of solving the ‘adjoint’ linear system, A T λ = c for λ T = c T A 1. Note that the original problem is assumed to be represented by the system of linear equations, Au = b. Hence, the gradient for each step of an optimization algorithm (we use quasi-Newton) requires a single extra linear solve, rather than a linear solve for each of the many parameters, p. This saving makes it possible to do large inverse problems, such as computing a basal traction for each point in the model domain (Fig. 2). Figure 3 corresponds to the initialized velocity and temperature field, or the steady-state solutions to the field equations that assimilate the data. This will provide the starting point for all numerical experiments. In some cases, such as determination of the sensitivity to q g, the entire assimilation/initialization process is repeated with different values.

Fig. 2. Top: modeled and observed velocity, as well as the portion of the modeled velocity accounted for by internal deformation. Middle: the β 2 field derived from the data-assimilation procedure. Bottom: the topography underlying the modeled ice profile.

Fig. 3. (a) Velocity and (b) temperature fields produced by the data-assimilation process. White lines in (a) indicate flowlines within the velocity field.

Numerical Experiments

Sensitivity of the FMB to basal heat flow

In order to determine the sensitivity of the location of the FMB to different values of geothermal heat flux, we performed the data-assimilation procedure over a range of possible values. This experiment is motivated by the observation that basal sliding represents a significant portion of the total modeled surface velocity, and we wish to determine the geothermal heat flux required to produce a completely melted bed, in line with the assumption that the bed must be at the melting point for sliding to occur. We conducted model runs every 5 mWm–2 within the range 0–120mWm–2. Figure 4 shows the location of the FMB as a function of the prescribed geothermal heat flux. The FMB asymptotically approaches the divide as geothermal heat flux is increased, although the entire bed is not at the melting point under any of the parameter values considered, even for fluxes which seem unreasonably high. For comparison, previous authors have used a value of 42 mWm–2 (Pattyn, 2003), and a structural similarity model by Shapiro and Ritzwoller (2004) indicates a mean geothermal heat flux of ~58mWm–2 along our flowline.

Fig. 4. Sensitivity of FMB location to variations in the geothermal heat flux.

Sensitivity of the FMB to sliding

Previous work suggests that seasonal changes in the glacial drainage system below the equilibrium-line altitude (ELA) can contribute to changes in basal traction, leading to changes in surface velocity (Zwally and others, 2002; Joughin and others, 2008; Van de Wal and others, 2008; Bartholomew and others, 2010). There is little agreement between these authors regarding the magnitude of proposed changes in surface velocity. Joughin and others (2008) suggest that terrestrially terminating glaciers in the region south of Jakobshavn Isbræ (of which Isunnguata Sermia is one) experience 25% increases in surface velocity as a result of surface meltwater lubricating the bed. Bartholomew and others (2010) suggest speed-ups as great as 200%, and that a warming climate and associated surface lowering will expose greater portions of the bed to surface meltwater, increasing the fraction of the ice sheet exposed to summer speed-ups. Van de Wal and others (2008) acknowledge these seasonal variations, but present data showing an overall 10%decrease in surface velocities between 1990 and 2007. They also note that surface ablation and velocity show no correlation.

Regardless of the magnitude and sign of such changes in surface velocity, we sought to determine whether perturbations to the basal traction field, β 2, downstreamfrom the ELA, such as those which would be induced by increased surface meltwater production, would have an impact on the basal thermal regime, specifically the location of the frozen melted boundary. We tested this by inflicting constant multiplicative perturbations to β 2 downstream from the ELA, ranging between 50% and 200% of the value produced by the data-assimilation process. The location of the FMB was insensitive to all these perturbations. The reason for this is shown in Figure 5. Notable changes in the surface velocity field (>1ma–1) extend only 20 km, or ~10 ice thicknesses, upstream from the extent of the perturbation.

Fig. 5. Sensitivity of surface velocity to perturbations to the basal traction field, β 2. ELA is ~1500ma.s.l.We interpret the longitudinal coupling threshold (LCT) to be the location at which the difference between any two surface velocity profiles is ≤1ma–1.

Thus, the advection of heat away from the bed, the dominant mechanism accounting for heat flux at the bed, as shown in Figure 6, is unchanged 90 km upstream, at the location of the FMB. This short coupling distance within the velocity field is corroborated by other studies (Price and others, 2008; Bartholomew and others, 2010).

Fig. 6. Budget of (a) heat sources and (b) sinks along the profile basal boundary. Latent heat generation (not shown) is a negative nonzero term below the FMB, and accommodates excess heat generated from (a). Strain heat is a positive nonzero term, but negligible compared to frictional and geothermal heat sources.

Heat budget

In order to track the dominant factors which dictate the thermal regime at the bed, we calculated a heat budget of sources and sinks in terms of flux to the ice-sheet base. We performed this calculation for a model scenario with optimized β 2 and a geothermal heat flux of 42 mWm–2; results are displayed in Figure 6. Upstream of the FMB, frozen conditions are controlled by the advection of cold ice.

Near the divide this advection is predominantly vertically directed from the surface. Moving downstream, advection becomes bed-parallel, so the advective flux decreases to zero at the FMB as heat sources raise the ice temperature as it flows along the bed. Throughout the frozen zone, and some kilometers beyond the FMB, the primary source of heat along the bed is geothermal heat flux. We find that heat generation due to straining at the bed is a positive contributor, but negligible compared to geothermal and frictional sources. Downstream of the FMB, excess heat generation is accommodated by the consumption of latent heat associated with the phase transition from ice to water as basal melt occurs. Basal melting initiates at the FMB and steadily increases to a maximum of ~20mma–1 near the terminus.


The sensitivity experiments described above indicate strongly different responses by the FMB to perturbations in geothermal heat flux and basal sliding. The direct response of the basal thermal regime to changes in geothermal heat flux is an expected result. However, the diminishing sensitivity of the FMB to increasingly higher heat fluxes is worth noting, and likely reflects the inability of the added heat to overcome cold advected from upstream.

In contrast, longitudinal coupling effects from sliding perturbations below the ELA do not propagate far enough up-glacier to influence the FMB. The location of the FMB is consequently insensitive to such perturbations. This interpretation hinges on the assumption that sliding perturbations apply only below the ELA, which is reasonable considering that the effect of increased surface melt input to the basal hydrologic system is not likely to propagate a significant distance upstream along the bed. We hypothesize that the limited distance over which longitudinal coupling occurs is a result of stress being dissipated at the basal boundary. It is important to note that sliding in our model is not limited to below the ELA. In fact, our optimization scheme produces a β 2 field with sliding above the FMB to the ice divide, albeit the upstream sliding is small relative to that occurring near the margin. Migration of melted conditions to the divide does not occur under very high values of geothermal heat flux, so the representation of sliding in our modeled frozen zone needs some explanation.

If the bed is in fact frozen we see several potential explanations for our modeled sliding. Firstly, sliding has been observed over a frozen bed consisting of a till layer (Engelhardt and Kamb, 1998), or hard bedrock (Echelmeyer and Wang, 1987; Cuffey and others, 1999). Additionally, substantial deformation has been observed within a frozen till layer, both within the body of the till itself (Echelmeyer and Wang, 1987; Engelhardt and Kamb, 1998) and along discrete shear planes (Echelmeyer and Wang, 1987). This mechanism may be taking place if such a layer exists beneath GrIS. Secondly, and perhaps more likely, our model could under-represent velocity from ice deformation, requiring our optimization scheme to over-represent sliding to maintain the observed surface velocity. Changes in flow due to variable impurity and water content and grain size of ice are not accounted for in our model; however, elsewhere in Greenland, a layer of soft pre-Holocene ice has been observed to enhance flows 1.7–3.5-fold (Paterson, 1994; Lüthi and others, 2002). Alternatively, but in the same vein, the standard constitutive law we use could under-represent grain-scale ice deformation processes (Goldsby and Kohlstedt, 2001). Finally, the velocity field itself could potentially depict a velocity field out of balance with the current ice-sheet geometry. We have no basis to eliminate any of these possibilities. However, if the magnitude of sliding over frozen bed computed here is not real, it is likely to be principally accounted for by spatial changes in geothermal heat flux, anisotropies within the ice, or a combination of the two.

An alternative scenario is therefore a partitioning of the observed surface velocity with enhanced ice deformation and reduced sliding velocity. Our heat budget along the basal boundary suggests the implementation of sliding has a significant influence on the location of the FMB by increasing the advection of cold ice along the bed. Under this alternative scenario, the associated drop in sliding velocity combined with additional interior heat generation from enhanced straining may modulate the cold contribution from advection, pushing the FMB further upstream. These processes would likely be countered by a decrease in frictional heating, which would force the FMB towards the margin. A model exploration of FMB migration from this interaction of heat sources is beyond the scope of this paper, but will be the focus of future work.


We have developed a thermomechanically coupled, two-dimensional flowline model and applied it to a terrestrially terminating glacier profile located in western Greenland. We extracted geometric information for the model domain from a dataset presented by Bamber and others (2001), which we believe represents the best data available at present for the Greenland ice sheet. We used adjoint methods to optimize the basal traction field, such that modeled surface velocities matched observed values (Joughin and others, 2010) to within 1 ma–1.

With an optimized model in hand, we conducted experiments in order to determine the sensitivity of the FMB to perturbations in the basal heat flux and basal traction downstream of the ELA. We found that the FMB migrates easily under different assumptions about geothermal heat flux. At values close to 0mWm–2, the FMB moves very close to the terminus, but part of the bed remains unfrozen due to frictional heating from sliding. At high values, the FMB asymptotically approaches the ice divide. We found that for reasonable perturbations to basal traction downstream from the ELA, such as might be expected from an increase in surface meltwater production and associated bed lubrication, the FMB was insensitive. This is a result of the short length scale over which longitudinal stress coupling in the ice operates (~10 ice thicknesses). The FMB is significantly further upstream from the ELA than the perturbations to the velocity field extend, so advective heat fluxes are unchanged. Our model predicts that, under most assumptions about geothermal heat flux, sliding occurs over a frozen bed. We see two possible explanations for this: (1) that this sliding is real, and follows one of the mechanisms proposed by Echelmeyer and Wang (1987), Engelhardt and Kamb (1998) or Cuffey and others (1999); (2) anisotropies or variability in hardness within the ice result in a preferential flow direction and increased deformation. Additional work is needed to quantify the sensitivity of the basal thermal regime to the second of these factors.


Funding for this work was provided by the Greenland Analogue Project (GAP), a collaborative project funded by the nuclear waste management organizations in Sweden (Svensk Krnbrnslehantering AB), Finland (Posiva Oy) and Canada (NWMO). It was also funded by grants from the US National Science Foundation: NSF-CMG-0934404 (J.V.J.) and NSF-OPP-0909495 (J.T.H.).


Andersen, K.K. and 50 others. 2004. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature, 431(1005), 147–151.
Bamber, J.L., Layberry, R.L. and Gogineni, S.P.. 2001. A new ice thickness and bed data set for the Greenland ice sheet. 1. Measurement, data reduction, and errors. J. Geophys. Res., 106(D24), 33,773–33,780.
Bartholomew, I., Nienow, P., Mair, D., Hubbard, A., King, M.A. and Sole, A.. 2010. Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier. Nature Geosci., 3(6), 408–411.
Cuffey, K.M., Conway, H., Hallet, B., Gades, A.M. and Raymond, C.F.. 1999. Interfacial water in polar glaciers and glacier sliding at –17◦C. Geophys. Res. Lett., 26(6), 751–754.
Dahl-Jensen, D. and 6 others. 1998. Past temperatures directly from the Greenland ice sheet. Science, 282(5387), 268–271.
Davis, T. 2004. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. AMS Trans. Math. Softw., 30(2), 165–195.
Deuflhard, P. 1974. A modified Newton method for the solution of ill-conditioned systems of nonlinear equations with application to multiple shooting. Num. Math., 22(4), 289–315.
Echelmeyer, K. and Wang, Z.. 1987. Direct observation of basal sliding and deformation of basal drift at sub-freezing temperatures. J. Glaciol., 33(113), 83–98.
Engelhardt, H. and Kamb, B.. 1998. Basal sliding of Ice Stream B, West Antarctica. J. Glaciol., 44(147), 223–230.
Ettema, J. and 6 others. 2009. Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modelling. Geophys. Res. Lett., 36(L12), L12501. (10.1029/2009GL038110.)
Fahnestock, M., Abdalati, W., Joughin, I., Brozena, J. and Gogineni, P.. 2001. High geothermal heat flow, basal melt, and the origin of rapid ice flow in central Greenland. Science, 294(5550), 2338–2342.
Goldsby, D.L. and Kohlstedt, D.L.. 2001. Superplastic deformation of ice: experimental observations. J. Geophys. Res., 106(B6), 11,017–11,030.
Greve, R. 2005. Relation of measured basal temperatures and the spatial distribution of the geothermal heat flux for the Greenland ice sheet. Ann. Glaciol., 42, 424–432.
Greve, R. and Hutter, K.. 1995. Polythermal three-dimensional modelling of the Greenland ice sheet with varied geothermal heat flux. Ann. Glaciol., 21, 8–12.
Hughes, J.R. 2000. The finite element method: linear static and dynamic finite element analysis. Mineola, NY, Dover Publications.
Huybrechts, P., Payne, T. and the EISMINT Intercomparison Group. 1996. The EISMINT benchmarks for testing ice-sheet models. Ann. Glaciol., 23, 1–12.
Joughin, I., Das, S.B., King, M.A., Smith, B.E., Howat, I.M. and Moon, T.. 2008. Seasonal speedup along the western flank of the Greenland Ice Sheet. Science, 320(5877), 781–783.
Joughin, I., Smith, B.E., Howat, I.M., Scambos, T. and Moon, T.. 2010. Greenland flow variability from ice-sheet-wide velocity mapping. J. Glaciol., 56(197), 415–430.
Lüthi, M., Funk, M., Iken, A., Gogineni, S. and Truffer, M.. 2002. Mechanisms of fast flow in Jakobshavn Isbræ, West Greenland: Part III. Measurements of ice deformation, temperature and cross-borehole conductivity in boreholes to the bedrock. J. Glaciol., 48(162), 369–385.
MacAyeal, D.R. 1993. A tutorial on the use of control methods in ice-sheet modeling. J. Glaciol., 39(131), 91–98.
Oswald, G.K.A. and Gogineni, S.P.. 2008. Recovery of subglacial water extent from Greenland radar survey data. J. Glaciol., 54(184), 94–106.
Paterson, W.S.B. 1994. The physics of glaciers. Third edition. Oxford, etc., Elsevier.
Paterson, W.S.B. and Budd, W.F.. 1982. Flow parameters for ice sheet modelling. Cold Reg. Sci. Technol., 6(2), 175–177.
Pattyn, F. 2003. A new three-dimensional higher-order thermomechanical ice-sheet model: basic sensitivity, ice stream development, and ice flow across subglacial lakes. J. Geophys. Res., 108(B8), 2382.
Price, S.F., Payne, A.J., Catania, G.A. and Neumann, T.A.. 2008. Seasonal acceleration of inland ice via longitudinal coupling to marginal ice. J. Glaciol., 54(185), 213–219.
Shapiro, N.M. and Ritzwoller, M.H.. 2004. Inferring surface heat flux distribution guided by a global seismic model: particular application to Antarctica. Earth Planet. Sci. Lett., 233(1–2), 213–224.
Strang, G. 2007. Computational science and engineering. Wellesey, MA, Wellesey-Cambridge Press.
Van de Wal, R.S.W. and 6 others. 2008. Large and rapid melt-induced velocity changes in the ablation zone of the Greenland Ice Sheet. Science, 321(5885), 111–113.
Weertman, J. 1968. Comparison between measured and theoretical temperature profiles of the Camp Century, Greenland, borehole. CRREL Res. Rep. 246.
Zwally, H.J., Abdalati, W., Herring, T., Larson, K., Saba, J. and Steffen, K.. 2002. Surface melt-induced acceleration of Greenland ice-sheet flow. Science, 297(5579), 218–222.