Introduction
The interactions of ice sheets with climate include processes which operate globally, on long time-scales, and for which few data are available. These interactions must therefore be studied with models. Given the wide concern about the possible warming effects of increased concentrations of gases that absorb infrared radiation at terrestrial frequencies, especially carbon dioxide, methane, and chlorofluorocarbons, we cannot ignore climatic forcing on glaciers. Climate-model simulations generally agree that Arctic temperature increases will be large, in the range of 6° to 14° C, but these large changes generally result from sea-ice changes and are concentrated over the winter oceans. Nevertheless, using the study of Reference Washington and MeehlWashington and Meehl (1984) as an example, the predicted temperature increases for the Canadian Archipelago are in the range of 4° to 6° C in winter and 2° C in summer. Models are in less agreement over changes in precipitation amounts. These changes are sufficient to have a major affect on Arctic glaciers. It is therefore of interest to develop means by which the effects of changed climate on glaciers may be evaluated. This paper uses a finite-element model of a flow plane of Barnes Ice Cap, Baffin Island, to evaluate the current energy balance of the ice cap and to investigate the sensitivity of the ice cap to surface-temperature changes.
Description Of The Finite-Element Model
The current model was developed by coupling and enhancing two earlier flow-plane models discussed by Reference Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979,1980). One of these models calculated steady-state temperatures, and the other calculated steady-state velocities. Both models held density, heat capacity, and thermal conductivity constant over the flow plane. Although these models used different types of finite-element procedure for their solution, each discretized the vertical flow plane into quadrilateral elements, so both could be used on the same flow-plane grids. In the current model, temperature and velocity calculations are done together using similar numerical procedures. Additionally, the temperature calculation has been made fully time-dependent, allowing the model to be stepped through time. Densities, heat capacities, and thermal conductivities, which were formerly held constant, now vary.
Part of the advantage of finite-element methods for glaciological applications comes from their lack of dependence on rectangular grid systems. This makes them more suitable for odd shapes and saves computation time by requiring high densities of nodal points only in areas where the independent variables are likely to have high gradients. A more subtle advantage of the method is that computer programs created for finite-element solutions are usually independent of any particular grid. The method of finding integral constraints, the geometry of the finite elements, and the assumptions used to approximate the variation of a field over the elements can all have significant effects on the feasibility of a calculation and even on the results produced. Sufficient information is included in the description of the present model so that comparisons between this model and others can take these differences into account.
The basic equations of the model are the usual ones for steady conservation of momentum, conservation of mass for an incompressible medium, and conservation of thermal energy, expressed here as
where τij
is the deviatoric stress tensor, p is pressure, gi
is the gravity vector, ui
is velocity, k is thermal conductivity, θ is temperature, C is heat capacity, 
is the strain-rate tensor, Q is the strain heating rate, t is time, and xi
is the position vector. The convention of summing repeated indices over their range is used here. Since the model is for a flow plane, the gravity vector and partial derivatives have no component in the x
3 direction with one exception: the compressive or extensive transverse strain-rate, ∂u
3/∂x
3, can be specified as non-zero, in which case it will be included in Equation (2) and in the calculations of viscosity and heat generation below. For all other calculations, the model is two-dimensional, with x
1 as a horizontal coordinate and x
jpointing vertically upward.
The effective viscosity, μ, is defined such that
The viscosity is calculated as
where n is a constant usually taken to be 3, 
is the effective strain-rate, 
and B is a function of temperature and crystal structure, among other things. Equations (4) and (5) taken together are the familiar Reference GlenGlen (1955) power law for flow. The effects of temperature and crystal fabric on viscosity separate B into two variables:
where
E is a crystal-structure enhancement factor (Reference LileLile, 1978), θk is temperature in Kelvin, and B* is expressed in Pa year1/n . Equation (7) was developed by Reference HookeHooke (1981) from a review of many published laboratory and field measurements; it follows the form of an Arrhenius relation except near the melting point. Equation (7) was intended for use when n = 3, which is normally guessed in the absence of contrary data and which is used in this study. Specification of the viscosity relation also provides a simple form for the viscous dissipation term of Equation (3):
Densities and gravity must be presumed known by the model. If enhancement factors and transverse strain-rates are not assumed to be 1.0 and 0.0, respectively, they must also be specified to the model. Heat capacity and thermal conductivity are calculated as functions of temperature using data from Reference BuddBudd (1969) for ice near pure ice density:
where θ is in degrees Celsius, the units of к are J/(m deg year), and the units of C are J/(kg deg). Equations (4) through (10) provide relations for all of the quantities in Equations (1), (2), and (3), and allow us to solve for the independent temperatures, velocities, and pressures.
The numerical solution of this model begins with the discretization of the domain. Breaking the domain into M elements, each element must be chosen such that density, pressure, and viscosity may be assumed constant within the element, and the remaining quantities vary roughly linearly across the element. The independent. variables for this model will thus be the N temperatures θi, the 2N velocities 
and the M pressures pe.
Integral equations result directly from Equations (1), (2), and (3) using the method of weighted residuals. The elements used here are linear Taig quadrilaterals, each defined by its four corner nodes (Reference Irons and AhmadIrons and Ahmad, 1980). Each corner node has a corresponding shape function, defined in natural coordinates for the element. The natural coordinate system uses the two coordinates (ξ,n) and defines each element as a square whose sides are the ±1 coordinate lines. The local shape functions, relative to a single element, may then be defined as
where (ξi , ni ) represents the natural coordinates of node i.
The model uses two-point Gauss–Legendre quadrature in each dimension of the natural coordinates to evaluate integrals over the elements. The model produces non-linear equations which can be solved iteratively using successive substitution of an initial guess. The solution method is computationally slow but robust.
For the numerical calculations, the model separates into two parts: a dynamic model and a thermodynamic model. Each of these sub-models requires its own boundary conditions. Boundary conditions arise naturally from surface integrals created in the weighted-residuals procedure. Velocity model boundary conditions may be a velocity component or a force component. For a polar ice sheet, frozen at its bed, the boundary conditions are zero stresses at the upper surface, zero velocities at the base, and a mixed boundary condition of zero horizontal velocity and zero shear stress at an ice divide. Temperature or a temperature gradient are specified at each boundary node for the temperature model. In these experiments, climatic temperatures at the upper surface, basal temperature gradients to provide geothermal heating at the base, and zero temperature gradients at the ice divide are specified.
Basal temperature gradients are difficult to estimate and may change with time, as noted by Reference Budd, Young and Robin deBudd and Young (1983). In the present program, a bedrock-temperature model calculates basal temperature gradients from one-dimensional finite-element models of temperature diffusion in the rock beneath the ice sheet. These vertical bedrock-temperature models are run beneath the basal nodes of the flow plane, giving a two-dimensional bedrock-temperature field, although horizontal conduction in the rock is ignored. The bedrock-temperature models use basal ice temperatures as upper boundary conditions and a constant temperature gradient as lower boundary condition.
The independence of a finite-element model on any particular grid enables one to test the model in idealized geometries against analytical solutions. Input decks for bench-mark testing of this model were developed based on Reference NyeNye’s (1957) velocity solution for a slab of ice on a slope, Reference Robin deRobin’s (1955) steady-state temperature solution for a column, and the transient temperature solution of Reference Hanson and DickinsonHanson and Dickinson (1987). Additionally, a standard heat-conduction solution from Reference Carslaw and JaegerCarslaw and Jaeger (1959) was used to test the bedrock-temperature model. In all of the model-verification studies, one can get at least three-digit agreement between theory and model if one carefully removes from the model any effects not included in the theoretical solutions (such as strain heating) and chooses elements that reduce discretization error.
This model and the models used by Reference Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979), Reference RaymondRaymond (1983), and Reference NixonNixon and others (1985) all use the same quadrilateral element formulation. The velocity model used here, however, differs from these predecessors in the type of linear shape function used and the method of integration over the element. Of these preceding studies, only Reference Hooke, Raymond, Hotchkiss and GustafsonHooke and others (1979) used a temperature model, and then only in the steady-state form. The shape functions and integration method used in that temperature model have been adapted to the entire calculation here. A much different approach to the velocity calculations has been described by Reference HodgeHodge (1985). Hodge used quadratic shape functions on elements defined by nine nodes in a three-by-three array. The primary advantage of quadratic elements over linear elements is that the equation storage space and solution time of the model are reduced. If computer storage and time are not serious constraints, linear elements are more flexible with particular regard to the number of nodes that must be placed in a vertical string and the manner in which varying characteristics of ice, such as density, can be isolated. The greatest need for flexibility in element size and placement is generated by non-homogeneous ice.
The Present Temperature Conditions Of Barnes Ice Cap
Barnes Ice Cap is a small ice cap situated on a plateau on central Baffin Island (Fig. 1). The 150 km long ice cap has two domes: a 62 km wide north dome and a 22 km wide, nearly circular, south dome. Margins of the ice cap are generally at an altitude of 500 m. So far as is known, there is little relief in the topography beneath the ice cap (Reference JonesJones, 1972). Ice thicknesses are generally less than 500 m. Although the British Baffin Island Expedition of 1950 performed initial surveys of the ice cap, intensive study of Barnes Ice Cap dynamics began in 1970 when a trilateration net was established on the south dome (Reference JonesJones, 1972; Reference HoldsworthHoldsworth, 1975). This net consisted of a series of surface stakes spaced every 500–1000 m along the 10 km flow plane, with stakes also placed transvere to the flow plane to form a series of overlapping square strain nets. Repeated surveys and partial surveys along this flow plane continued until the spring of 1984. Hooke also drilled a series of short bore holes for near-surface temperature measurements, and a series of deep bore holes for measurement of the variation of temperature and strain-rates with depth. Data from the flow plane have supported previous work on its dynamics, thermodynamics, and crystal fabrics, including work with the previous versions of the present finite-element models (e.g. Reference Hooke, Raymond, Hotchkiss and GustafsonHooke and others, 1979,1980,1987; Reference Hooke and HudlestonHooke and Hudleston, 1980; Reference Hooke and HansonHooke and Hanson, 1986).
Map of a part of the South Dome of Barnes Ice Cap showing surface elevations above sea-level and positions of the bore holes and flow plane used for this study.
Interpolations of the surface and bore-hole temperature measurements approximate the current temperature field of the flow plane (Fig. 2). The time derivatives calculated from the finite-element model may then be used to evaluate the current state of energy balance on the flow plane. Additionally, the model may be used to calculate a steady-state temperature field which matches the current boundary conditions. Both of these calculations require a velocity field that is either input to the model or calculated by the model. The velocity field used here has been calculated by the finite-element model on a grid with 693 nodes and 588 elements after tuning to match the complicated crystal-structure variations on the flow plane.
Difficulties in understanding viscosity variations in the flow plane were dicusssed by Reference Hooke and HansonHooke and Hanson (1986). Two main types of ice comprise the flow plane: soft, bubbly white, Pleistocene ice primarily in a thin layer near the base, and less bubbly blue ice of Holocene age
Current flow-plane temperatures for the Barnes Ice Cap flow plane. Isotherms at 0.5° intervals, constructed from surface measurements and bore-hole measurements made during 1976–79. described bv Reference Hooke and HudlestonHooke and others (1980).
throughout the rest of the flow plane. Additionally, the blue ice has several layers of varying crystal fabrics. To successfully simulate the known velocity field at the surface and in the bore holes, it has been necessary to segregate the white ice and assign it an enhancement factor, E, of 5.8. Additionally, a layer of blue ice exhibiting multiple maxima of preferred crystal orientation has been stiffened with E = 0.6. Strain-rates in the rest of the ice cap are successfully simulated with E = 1.0.
Basal temperature gradients at each bore hole were calculated by Reference Hooke and HudlestonHooke and others (1980) for use in a one-dimensional finite-difference model of temperatures in the bore holes. These temperature gradients were partly based on the observed gradients in the measured temperature profiles near the bottom of the bore holes, and partly on theoretical considerations of obtaining a nearly steady-state solution for the model. Their pattern indicates a response in the bedrock to warming temperatures at the base, especially near the margin. These data have been interpolated over the flow-plane base to provide initial conditions to the rock-temperature model. They have also been interpolated with depth to match the inferred deep rock-temperature gradient of −0.015° C m1, according to the formula
where y is the vertical coordinate of a rock column of thickness Η, β is the vertical temperature gradient, and the base (deep rock) and top (ice-rock interface) of the bedrock column are denoted by the subscripts B and T, respectively. Equation (12) is the analytical solution for diffusion within the rock with temperature gradients specified at both ends of the column, using only the steady-state term and the first eigenfunction.
The current energy balance of the flow plane can be understood somewhat by comparing averages of velocity components and temperatures over the flow plane under current conditions and under several different steady-state assumptions (Table I). All of these steady-state calculations use the current surface-temperature boundary conditions. The first two compare the effects of using basal temperature gradients calculated by Reference Hooke and HudlestonHooke and others (1980) and using a constant basal temperature gradient of −0.015°C m−1. The latter case gives a higher input into the base of the ice cap, leading to higher temperatures and especially higher velocities, as the greatest temperature increases are in the region of highest strain near the bed.
The last two steady-state calculations in Table I show the effects of calculating a temperature field that is in steady state with the current temperature field and the two choices of basal temperature gradient, and is also in steady state with the current velocity field. These are not realistic steady states, but they show the mutual influence of temperature and velocity. Allowing velocity increases produces a slightly cooler flow plane, since the average direction of velocity is from the coolest parts of the ice cap. These effects are very small, such that one could say that velocity is more sensitive to temperature than temperature is sensitive to velocity.
Flow-plane average temperatures and velocities for steady-state conditions on the present barnes ice cap tril ater ation net flow plane. “Temperature steady state” indicates that velocities were held constant, and βT is the temperature-gradient boundary condition at the base of the ice
To illustrate the approach to steady state, the model was run with current initial conditions for 2000 years of simulated time in 10 year increments, with the bedrock-temperature calculation included (i.e. the model should approach the first steady-state condition of Table I). The temperature curve overshoots steady state during the run because of different temperature trends in different parts of the flow plane (Fig. 3). The initial temperature field contains a cooling trend near the divide and a warning trend throughout most of the flow plane. In the short run, these conditions lead to a net warming over the flow plane.
Average flow-plane temperatures for the Barnes Ice Cap flow plane in a 2000 year model simulation using the present temperature field as initial condition and holding surface-temperature boundary conditions constant at current values.
In the long run, the cooling near the divide decreases the horizontal temperature gradients through the central area of the flow plane, so that the cooling area is effectively advected into a larger area of the flow plane. Since this expansion of the cooling area occurs as the temperatures in the thinner ice near the margin are beginning to approach final steady-state values, the temperature of the flow plane begins to decrease.
Sensitivity To Forced Temperature Changes
A primary purpose of this study has been to understand how surface-temperature changes affect the internal conditions of an ice cap and the response times to these conditions. Barnes Ice Cap has served as a starting point for these calculations because of data availability and because it may be considered typical of small Arctic ice caps. For the following calculations, an input deck with idealized geometry and simplified temperature structure was produced in order that the results would be less influenced by the particular characteristics of Barnes Ice Cap. The base was idealized to a horizontal surface, and the upper surface is a fifth-order polynomial which was fitted to the actual surface profile by least squares. The upper-surface profile is visually indistinguishable from the original when plotted, but the idealized version has a substantially smoother slope which increases monotonically towards the margin. All enhancement factors were set to one, and all transverse strain-rates were set to zero. The grid consists of an array of elements with eight in the vertical and 249 in the horizontal. The nodes are evenly spaced in the horizontal direction. The vertical thicknesses of the elements increase towards the interior from 15 m thick elements at the surface and the base.
The basal temperature gradient was set to −0.015°C m−1, and the surface temperature profile was set to follow the U.S. Standard Atmosphere lapse rate of 6.5°C km−1. Stream lines and isotherms of the resultant steady state are shown in Figure 4. Flattening the base of the ice cap produced greater ice thicknesses through most of the flow plane, which would lead to greater velocity except for the removal of the enhanced white ice at the base. The resulting average velocity components are 3.09 m year−1 in the horizontal and −0.120 m year−1 in the vertical. The flow-plane temperature has been cooled to —10.5°C by the assumed lapse rate. The stream lines produced are qualitatively similar to the original ice cap, with the equilibrium line moved slightly up-glacier. The distortions in the isotherms are as expected from temperature advection along the stream lines.
The magnitude of the response to thermal forcing was first investigated by increasing the temperature evenly over the surface and calculating the steady-state response of the flow plane to each higher surface temperature. These calculations were carried out for temperature increases up to 6°C. The response of flow-plane temperature is linear but not direct (Fig. 5). For each degree of surface-temperature increase, the average flow-plane temperature increases 0.9°C, and this ratio does not change within the limits of accuracy of these calculations. Velocities increase in response to rising temperature, and their increases are non-linear, as would be expected from Equation (7). As noted earlier, velocity changes for a given temperature response are larger relative to their intitial sizes than the corresponding temperature changes.
There are two feed-back mechanisms by which velocity changes prevent the flow-plane temperature from directly following the surface temperature. One is a positive feedback, as increased velocity leads to increased strain-rates, increased strain heating, and further increases in temperature. The second is negative, as increased velocities lead to additional advection of ice from the ice divide and the upper surface, which are the coldest parts of the ice cap. Since the response of the ice cap is smaller than the temperature forcing, the second feed-back is more important. The relative sizes of these effects were estimated from steady-state calculations in which strain heating was set zero. The flow-plane temperature with current boundary conditions dropped to −10.81°C, and the change in response to a 1°C forcing dropped to 0.886. The positive feed-back of strain heating is thus about 12% as large as the negative feed-back of advection.
Steady-state temperature and velocity solutions for the idealized Barnes Ice Cap flow plane, (a) Isotherms of the temperature solution; contour interval is 0.5°C. (b) Stream lines of the velocity solution.
The net feed-back for the entire flow plane is —10%, and the feed-back at individual points is in the range of −0.8 to −0.17 over most of the flow plane below 80 m in depth and further than 2 km from the margin. The surface has a feed-back of zero, forced by boundary conditions, and a net positive feed-back as high as 28% occurs near the margin in the area that has an upward velocity component. The positive feed-back area mostly results from the high strain-rates near the margin. When strain heating is suppressed, positive feed-backs are confined to a small area within 1 km of the margin with a maximum of 8%.
The layer of soft, white ice near the base in Barnes Ice Cap is believed to be common in Arctic ice caps. The possibility that a more deformable basal layer changed these feed-backs was investigated with an additional model calculation. The basal layer of elements in the idealized deck is 15 m thick from the ice divide to 9370 m from the divide. At that point, the total ice thickness becomes less than 120 m, and all element thicknesses become one-eighth the total thickness. Enhancement factors in the bottom element layer were increased from 1.0 to 6.0 in increments of 0.1. At each enhancement-factor level, the flow-plane average temperatures and velocity components and their responses to 1 C temperature increases were calculated. As the enhancement factor increased from 1.0 to 6.0, the average velocity components slightly more than doubled in size (3.1 to 6.3 m year−1 for horizontal velocity, −0.12 to −0.26 m year−1 for vertical velocity). The changes in velocity were almost exactly linear with the change in enhancement factor.
Responses to the forced 1°C surface warming decreased from 0.900°C to a minimum of 0.892°C for E = 3.0, followed by an increase to 0.896°C. The effect of removing strain heating is to decrease the response by 0.012° to 0.014° C, with the minimum effect also occurring at E = 3.0, but followed by a constant trend rather than an increasing trend. The dominant change in feed-back comes from the negative feed-back effects of advection, which increase as the velocity increases, in accordance with expectation, for moderate basal softening. The change in trend for larger basal enhancements implies that as the strain becomes more concentrated at the base, the stream lines become more coincident with the isotherms. Because of the smallness of the changes in feed-back, it is difficult to demonstrate that this is in fact what is going on.
Changes averaged over the idealized flow plane as the surface-temperature conditions boundary are uniformly increased from 0° to 6°C expressed as departures from the 0°C case, (a) Temperature; (b) Horizontal component of velocity; (c) Vertical component of velocity.
To determine the response of the flow plane going from one steady-state condition to another, four experiments were run with the initial condition of the idealized grid subjected to an instantaneous temperature rise at the surface. These cases are the combinations of two temperatures changes, 1° and 2°C, with the inclusion of bedrock temperature-profile changes (Fig. 6). The bedrock slows the response proportionally to the size of the perturbation, such that at the end of 2000 years of simulated time, the ice-only simulations are within 10% of their final values,
Transient response of temperature, averaged over the idealized flow plane, in which the initial state shown in Figure 3 responds to a uniformly applied surface-temperature increase. Dashed curves show the response when bedrock-temperature variations are suppressed; solid curves show the response of the complete model. Responses for 1° and 2° C temperature forcing are shown.
while the simulations that include bedrock-temperature changes are more than twice as far from steady state. The faster response near the margin is illustrated by isochrones of the basal temperature change (Fig. 7).
The effect of including bedrock in the temperature calculations is more apparent if one calculates the response times implied by Figure 6, assuming an e-folding response time τ between each time step in which
where Δθ are the departures of temperature from steady state. When such calculations are applied to a complex system, the response time will be an average reponse time over all the eigenfunctions of the solution. Short response times correspond to high eigenvalues and small-scale spatial variations, and these eigenfunctions become unimportant rapidly because of their short response times. The average response time of the flow plane should thus increase asymptotically towards a maximum that corresponds to the lowest eigenf unction of the system. That fairly well describes the behavior of the model solution when the bedrock-temperature calculation is not included, as the response time reaches a steady value of 830 years within
Isochrones of basal temperature change for a 1°C increase in surface temperature on the idealized Barnes Ice Cap flow plane, at 500 year intervals labelled in centuries.
Response-time variations implied by the temperature variations in Figure 5. Solid curve is for the complete model; dashed curve is for suppressed hedrock-temperature calculations.
500 years of simulation (Fig. 8). Reference Hanson and DickinsonHanson and Dickinson (1987) determined from one-dimensional calculations that the maximum thermal-response time for the accumulation region is approximately given by
where H is the thickness and A is the accumulation rate. Applying this formula to the main accumulation zone of the idealized ice cap gives maximum response times mostly around 1100 years, leading to the conclusion that the two-dimensional response is faster than the one-dimensional response, but the one-dimensional estimate is correct to order of magnitude.
The second case in Figure 8 includes bedrock temperatures in the model calculation, but not in the average temperatures from which these response times were calculated. Early on, the curve is clearly dominated by changes within the ice cap. Forcing of the bedrock only begins when the temperature changes have begun to affect the basal temperatures of the ice, so increases in response time caused by the thermal inertia of the bedrock become important later. There is no clear evidence of a leveling off in this solution. The maximum response time of a one-dimensional column of rock under diffusive heat transfer is given by
where H is now the thickness of the rock column and it is its thermal diffusivity. For the current 4000 m thick column using a thermal diffusivity of 39 m2 year−1, we get a maximum response time of 42 000 years. This substantiates the conclusions of Reference Budd, Young and Robin deBudd and Young (1983) that one cannot ignore the bedrock-temperature profiles in long simulations.
Discussion And Conclusions
The model experiments described here examine two aspects of the temperature response of an ice cap to climatic thermal forcing: the timing of the response and the size of the response. Thermal response times of even this small ice cap are of the order of 1000 years, while significant climatic variations occur an order of magnitude faster (e.g. Medieval warm period followed by the Little Ice Age followed by the present warm period followed by enhanced greenhouse-effect warming). Discrepancies in response time produce such unexpected effects as the change in direction of response shown in Figure 3, indicating that Barnes Ice Cap is currently responding to more than one past climatic variation. It is worthwhile to remember when considering anticipated changes in Arctic glaciers that these glaciers may still be responding to changes associated with the Little Ice Age and thus exhibit a complicated response to additional forcing.
The transient responses shown in Figure 6 could be misinterpreted as indicating that bedrock-temperature calculations are relatively unimportant for such short-term temperature changes as may occur over the next century, as the effects of basal temperature-gradient variations do not manifest themselves for several centuries of simulation. The most important temperatures in an ice cap are those nearest the base which determine the viscosity in the highest strain region and which determine whether or not the glacier is frozen to its bed. This latter effect in particular will be important to ice caps in the marginal Arctic. Basal melting may ensue in some parts of Barnes Ice Cap during a significant warming, in which case the responses shown here would presumably become a small part of the effects. Because the basal temperature response is fastest near the margin, a receding margin would be more likely than a catastrophic onset of basal melting over the entire ice cap. The response of the basal temperatures will be fastest in the shallow ice near the margin, especially in the area where velocities have an upward component. The bedrock-temperature response is forced by basal temperature changes in the ice cap, so the response time for bedrock reflects the response time in a particularly interesting part of the ice cap.
The size of the temperature response approximates the size of the forcing, with modifications from two effects of the velocity response. The negative feed-back of increased advection is about 12%, while the positive feed-back of increased frictional heating is about 2%, leading to a net negative feed-back of 10%. Both feed-backs are closely tied to basal strain-rates, such that increasing basal strain does not change the relative importance of the feed-backs because it affects both similarly. Despite the negative feed-back of the flow plane as a whole, a small area near the base and margin exhibits feed-backs as high as 28%. The positive feed-backs occur in a small but significant region which has the highest temperatures and hence the greatest potential for the onset of basal melting.
Changes in mass balance have not been accounted for in this study. By running long simulations without change in the shape of the ice cap, we implicitly assume a steady balance between the emergence velocity at each point on the surface of the ice cap and the net mass balance at that point. I do not wish to assert that such a balance is likely to be maintained on average over the 2000 year periods shown in some of these simulations. Rather, I have used the computationally simplest approach to a situation in which no good prediction of mass balance is possible.
The difficulties with mass balance arise from the uncertainty of climate-model results. Changes predicted for the Arctic are small enough to be of uncertain sign. Inter-annual variations of net mass balance for Barnes Ice Cap are most strongly influenced by summer temperatures which are indicators of the length of the melt season, but that influence has been detected because of the relative smallness of inter-annual variations in precipitation (Reference Hanson and DickinsonHanson, 1987). The conventional wisdom that warmer winters have more snowfall results from higher frequencies of penetration of maritime air masses in some years, rather than on any direct, physical correlation between temperature and likelihood of precipitation. One could use both of these inter-annual trends to assert that, under global warming, Arctic glaciers will experience higher winter snowfall and higher summer melt. In this case, the sign of the change in net mass balance may be unknown, but the changes might be expected to be small because of the opposite effects.
The models described here are capable of simulating a transient response with prescribed mass-balance information and a continuously changing upper surface profile. In such experiments, velocity adjustments in response to changing flow-plane shape are much larger than in response to changing temperature fields, producing greater advective changes in the energy-balance calculations. It is impossible to focus on the dynamics of the ice and how it reached its current shape without taking such changes into account. Since these responses are on much smaller time-scales than the thermodynamic response and are of secondary importance in understanding the thermal response, they are beyond the scope of this study.
Popular accounts of global warming have in some cases drawn a direct correlation between rises in atmospheric CO2 concentrations and flooded shorelines, with the melting of a substantial part of the cryosphere as an intermediate step. The results in this paper cannot be extrapolated to temperate mountain glaciers nor to the two large ice sheets, but they should have some validity when applied to polar ice caps of the scale found on many of the Arctic islands. For these areas, the model-predicted temperature changes for the melt season are small, the time required for climatic temperature changes to influence greatly the basal temperatures and thus lead to possible melt and higher strain is of the order of centuries, and the overall change in ice temperature is slightly less than the surface-temperature change. While anthropogenic warming is likely to have a dramatic effect on the Arctic environment, its ultimate effects cannot be determined from simplistic reasoning, and the lag time will be long.
Acknowledgements
R.LeB. Hooke provided the Barnes Ice Cap data used in this study and made available the original versions of the finite-element models. The National Center for Atmospheric Research, which is sponsored by the U.S. National Science Foundation, provided computer support for the development of the models used here.
