Symbols
Introduction
A large deposit of iron ore has been discovered at Isua along the edge of the Greenland ice sheet (see Fig. 1). The feasibility study for the recovery of this ore involves many glaciological problems because the ore is mostly covered with ice up to 200 m deep. It is anticipated that the development of a large open-pit mine with ice walls around 60% of the perimeter will be necessary. As the ice overburden is stripped away to expose the ore, the flow of ice toward the open pit will increase as the surface slope is increased and the elevation of the ice margin is lowered. The volume of ice flowing toward the pit at the final stage of development is calculated as part of a preliminary cost analysis; the procedure and results are presented here.
Fig. 1. Map of Greenland showing location of Isua.
Many of the problems in this study are familiar to glaciologists but some innovations will be necessary to handle a project of this size. Also, some significant glaciological phenomena such as the occurrence of liquid water beneath the ice sheet will be investigated during the course of the feasibility study. A network of stakes for determining surface velocities has provided information for several years, but much additional information is needed in order to obtain surface velocities further up-stream from the ore. Likewise, the configuration of the glacier is known in the immediate vicinity of the deposit, but additional information must be obtained up-stream from the proposed excavation. Snow drifting along the margin of the glacier occurs at present and it is anticipated that the accumulation of snow in the open pit could be a severe problem. Many millions of cubic meters of ice will have to be removed directly, necessitating the development of rapid and economical methods of excavating and hauling ice. The excavation costs would be reduced if large amounts of ice could be removed by modification of the surface albedo with rock dust or carbon black. More information about the application and benefits of these materials must be obtained before the most efficient balance between direct excavation and surface melting can be found. The surface run-off resulting from melting will have to be controlled in order to prevent flooding of the mine, and large channels will be used to collect and remove water from the glacier surface. The stability of these water courses in ice over tens of years will be studied closely.
Although the mineral deposit has been sampled by drilling through the glacier, no ice samples were taken. Furthermore, temperature measurements in the ice taken to a depth of 115 m suggest that the melting temperature occurs around 200 m depth. A drilling program is planned to obtain ice samples for fabric analysis, to measure ice temperatures to greater depths, and to locate the subglacial water. Clearly these parameters have a significant effect on the rate of flow of the glacier, and, if a large area of basal water exists beneath the ice sheet as suggested by Gudmandsen and Reference Gudmandsen and ChristensenChristensen (1972), the mine could be flooded if the excavation intersects this subglacial reservoir.
If the results of the feasibility study are favorable, a new configuration will be imposed on the ice sheet over an area of about 3 km2. In a sense this will be the largest and most expensive glaciological experiment ever performed; glaciologists will have the opportunity to observe the behavior of a large ice mass which has undergone an artificial alteration. The planning of the excavation will require a series of compromises between the need for efficient mine development and the need to minimize the costs of ice excavation. Many years of research on ice mechanics and polar glaciers provide the basis on which the decisions will be made.
Physical Setting
The ore deposit at Isua is located along the margin of the ice cap as shown in Figure 2. The deposit is 2 km long and 250 m to 450 m wide with the long axis oriented approximately in the north-south direction. About 15% of the area of the ore is exposed at the surface including a rock ridge which rises some 60 m above the edge of the ice. The ridge causes flow diversion to the north and south of the ore body, thus flow rates are currently quite small in the center of the deposit (see Fig. 3) and the surface is nearly flat (see Fig. 4).
Fig. 2. Location of the ore deposit at the edge of the ice sheet. Regional velocities show the diversion of the flow away from the deposit.
Fig. 3. Surface velocities measured in the area of the ore body.
Fig. 4. Base map showing surface contour, bedrock contours and locations of the seven cross-sectional profiles. The outline of the ore deposit and the suggested areas of ore sacrifice are shown.
The regional flow (as determined by the surface contours) is from the east-north-east. Most of this flow is diverted to the north of the deposit by a subsurface ridge which trends east-south-east from the northern tip of the deposit. The flow rates to the north are as large as 30 m a−1 while those to the south are less than 15 m a−1. The flow rates at the center of the deposit are as small as 1 m a−1 where ice flows up-hill along the steeply sloping bedrock and is tapped by the ore deposit (see Fig. 5). The ice flows up-hill and over the deposit along its southern and northern flanks (see Figs. 6 and 7). The main outlet of the ice sheet in this area is the large discharge glacier which flows past the northern tip of the deposit and moves down-hill to the west Numerous smaller discharges occur to the south (such as the flow across the southern tip) but the large glacier on the north accounts for most of the ice flux out of this area.
Fig. 5. Cross-sectional view of the glacier along line No. 4 showing the bedrock, existing ice surface and proposed ice surface.
Up-hill flow occurs essentially from the center of the Greenland ice sheet (see Fig 8) with an average slope of less than 0.4% over the last 150 km. At Isua the slope increases to as much as 50% because of the presence of the resistant rock formations associated with the ore deposit The ridge of sub-surface rock on the north and the up-hill slope of the bedrock are favorable factors in reducing the potential flow into the open pit. No similar problem of up-hill flow has been analyzed, but the general principles of the theory of glacier flow can be applied to the problem treated here—the case of up-hill flow where the surface slope will be imposed on the glacier by direct excavation.
Fig. 6. Cross-sectional view of the glacier along line No. 6 showing the bedrock, existing ice surface and proposed ice surface.
Fig. 7. Cross-sectional view of the glacier along line No. 1 showing the bedrock, existing ice surface and proposed ice surface.
Fig. 8. Cross-sectional view of the Greenland ice sheet up-stream from Isua showing the inclination of the bedrock.
Flow Model
The complicated geometry of the glacier at Isua is very difficult to model and many simplifying assumptions are necessary. Accordingly, the area around the perimeter of the pit is divided into segments each of which is approximately by a two-dimensional glacier. The bedrock surface is taken as a sloping plane and a simple configuration is adopted for the proposed ice surface (see Fig. 9),
Although the adoption of a parabolic profile may preclude the selection of a configuration with lower flow rates, it is necessary to adopt a configuration which can be excavated easily. The parabolic profile, which is common in nature, may satisfy both criteria.
The development of the two-dimensional flow of a glacier begins with the equilibrium equations (Reference BuddBudd, 1969, p. 112),
Fig. 9. Idealized view of the existing and proposed glacier profiles.
and
The longitudinal variations must be retained in this problem, therefore the equations are treated by averaging the normal stress over the thickness of the glacier. Equation (2) is integrated over the thickness of the glacier and from the margin to x; σxz is taken as zero at the surface and σx is zero at the origin (see Fig. 9), so that
where σx is the average value of σx in the z direction. Likewise, Equation (3) is integrated from the surface (where σz is taken as zero) to z and then integrated over the thickness of the glacier,
where σz is the average value of σz in the z direction. For two-dimensional flow where no deformation occurs in the third direction, the averaged deviator stress in the longitudinal direction 
is related to the averaged normal stresses by (Reference BuddBudd, 1969, p. 214)
These last three equations are combined to give
In this equation the term on the left is equal to twice the force causing longitudinal deformations, the first term on the right describes the compression which arises due to up-hill flow, the second term describes the drag as the ice moves over the bedrock, the third term describes the longitudinal forces due to the curvature of the ice surface, and the final term describes other longitudinal variations in the shear stress.
A simplifying assumption is made in order to use Equation (7); the shear stress at the base is taken as a constant for all values of x. There is considerable justification for adopting this assumption since calculated values of basal shear stress are generally within 50% of 105 N m−2 (Reference KambKamb, 1964). The unusual geometry of flow at Isua suggests that an area of “dead material” could exist because the glacier is pushed against the bedrock. However, observations of the movement of drill holes suggest that the usual zone of deformation occurs at the bottom of the ice and that shear stresses are a maximum there. Using the assumption of a constant basal shear stress, Equation (7) is solved for σxz(x, o) at x = L and, upon substituting this expression for σxz (x, o) back into Equation (7),
The last two terms are dropped from this equation so that the stress can be calculated directly from any specified configuration. Budd (Reference Budd1969, p. 117) discusses the importance of these terms; their exclusion may cause some error in the calculation of the flow near the margin.
The flow law for ice is adopted from Reference GlenGlen (1955) in the form
where the values of B and n are chosen later according to the particular situation. The total thickness of the ice is approximately given by
Equations (8), (9) and (10) are combined to give
The validity of this equation is checked by examining the case of laminar flow 
over a level bed (β = o). The solution of Equation (11) gives m = 2, the classical profile for laminar flow.
The distribution of longitudinal strain-rate ėx can be calculated from this equation once the configuration has been specified and the value of the average deviator stress 
has been chosen at x = L. The distribution of velocity can then be calculated starting with a chosen value at the position L. In this way the ice flux at each point in the cross-section can be computed from values assigned to the glacier at one position. Ideally this position (L) would be located far enough away from the margin to ensure that the present values of velocity and stress would not be affected by the excavation. However, the practical limitations of obtaining the necessary information about the glacier over these distances compromise this principle. Therefore the equation of continuity,
where the total ice flux is
is used to select appropriate parameters at L. The temporal derivative disappears in the problem being considered here, the steady ice flux at the final stage of pit development, hence at the position L
The ablation rate ȧL at L is taken as 1 m a−1 net ice ablation, a value measured during only one field season. The thickness ZL and surface slope ∂Z/∂x|L are known for any suggested profile (a flat bed is assumed), therefore the longitudinal strain-rate and velocity are related through Equation (14). The use of this technique limits the choice of these parameters at the reference position, although much subjectivity still remains in their choice. For most of the cross-sections several combinations of these values were tried before a self-consistent result was obtained. This method of choosing the initial values is not entirely satisfactory, although little more can be done until more information about the geometry and flow up-stream from the ore deposit has been obtained.
Flow Calculations
The rock edge of the proposed pit is calculated from the outline of the ore body at the elevation of 885 m using pit walls sloping at 55 deg. The edge of the area of ice excavation (the reference position L) is chosen as the furthest position up-stream for which information is currently available. Generally this distance is 700 to 800 m from the edge of the pit, which, since we have the practical limitation of only being able to excavate economically over a limited distance, is a reasonable compromise among the different requirements placed on these locations. The flow toward the proposed pit is analyzed along the seven lines shown in Figure 4; five of these are expected to be principal directions of flow and two are lateral flow into the ends of the pit. Some flow from the eastern side of the pit will occur but it is negligible. The flows along three of the seven lines are discussed here and more details of this study are given in Reference GolbeckColbeck (1973).
A coordinate system is located at the proposed margin (see Fig. 9) and the bedrock surface is approximated as a sloping plane from the margin to the position L. The values of the parameters B and n in the flow law are chosen from Budd’s (Reference Budd1969, p. 23) compilation of laboratory and field evidence. The average temperature at the bedrock is taken at – 2° C from measurements in bore holes and the chosen values of n vary between two and four depending upon the overall level of effective stress. These choices are somewhat arbitrary in view of the controversy over the exponent in the flow law but they provide an upper limit on the rate of deformation. Clearly the flow would be much larger than that predicted here if liquid water exists at the base and additional flow occurs due to sliding. Current information suggests that the glacier is frozen to the bedrock over most of the area of interest and, since the excavation should take many years to complete, further freezing should occur at the rock–ice interface.
At each cross-section of flow, a trial value of m is chosen (this specifies the surface profile along that line) and the values of stress and velocity for the reference position are chosen subject to the constraint of Equation (14). The velocity, strain-rate and “surface flux” (the ablation plus excavation necessary to maintain the specified profile) are calculated for the entire profile. The value of m is then adjusted to eliminate undesirable patterns in the flow field and, once a satisfactory value of m has been found, the configuration of the ice surface and the total surface flux has been determined for that flow line. The total surface flux is the amount of ice which must be removed each year in order to keep the glacier from thickening and advancing into the mine. Much of this material can be removed by surface melting but continuous excavation will also be necessary.
The up-hill flow along line No. 1 is rather slow at present because of the subsurface rock ridge. The reference position L for the start of the flow calculation is set 500 m from the margin, the longitudinal velocity UL is chosen to be – 10 m a−1, and the longitudinal strain-rate ėL is chosen to be extending at 0.0075 a−1. The strain-rate increases rapidly toward the margin reaching a maximum value of 0.0138 a−1 at about 37 m (see Fig. 10). Reference MeierMeier (1958) found that crevasses form when extending strain-rates exceed 0.01 a−1, hence some crevassing is expected along the margin of this proposed profile. The surface flux equals the measured ablation rate (–1 m a−1) at the start of the profile but decreases very rapidly toward the terminus. Most of the excavation will be necessary in the region near the margin where high strain-rates, large surface slopes and some crevassing occurs. This may make the excavation difficult.
Fig. 10. The calculated values of strain-rate, velocity and surface flux (or ablation) along line No. 1.
The ore body reaches its lowest point between lines No. 4 and 5 and much flow will occur toward the proposed pit in this area. This study assumes full recovery of the mineral deposit, although it is recommended that two large areas of ore should be sacrificed in order to minimize ice movement. These areas are shown on Figure 4, the first being an area associated with the lowest elevation of the bedrock and the second being the northern area over which flow rates are currently quite large. The longitudinal strain-rate at the reference position L on line No. 4 is extending at 0.01 a−1 and the longitudinal velocity is –15 m a−1. Characteristically the strain-rate increases rapidly toward the terminus reaching a maximum value of more than 0.03 a−1 at 100 m from the margin (see Fig. 11). The surface flux decreases to 20 m a−1 at that position and then decreases more rapidly to –170 m a−1 at the margin.
Fig. 11. The calculated values of strain-rate, velocity and surface flux (or ablation) along line No. 4.
These high flow rates and surface fluxes occur because of the large decrease in the elevation of the terminus in this region (see Fig. 5), requiring the sacrifice of a large amount of ore as described above. In addition to the cost of the excavation, the dangers involved in working in an area with large extending strain-rates must be considered. The formation of large crevasses and seracs is likely on this steep glacier surface and avalanche hazards are high. It is possible that avalanches from collapsing seracs could reach the pit floor thus making working conditions very hazardous throughout this area. With the sacrifice of a small amount of the ore, the elevation of the terminus can be raised considerably thus reducing the problems associated with high flow rates.
At the southern end of the ore body the elevation of the pit wall is somewhat higher (see Fig. 6) thus reducing the flow toward the pit. The longitudinal velocity is taken as –6 m a−1 and the extending strain-rate as 0.003 a−1 at the reference position L. The strain-rate increases to more than 0.016 a−1 at 60 m from the margin and the surface flux decreases rapidly over the last 100 m (see Fig. 12). The longitudinal strain-rate is large enough to suggest that crevassing may occur but that it will probably not be a serious problem. If necessary the flow could be reduced further by leaving a small amount of ore in order to raise the margin of the glacier in this corner.
The details of the calculations and the results for the other four profiles are given by Reference GolbeckColbeck (1973). The total excavation necessary to establish the new surface profile and the total ice flux toward the open pit is calculated by summarizing the results from the seven flow lines The initial excavation requires the direct removal of 172 × 106 m3 of ice while the yearly ice flux into the excavated area is 7.9 × 106 m3 a−1. Much of the ice flux can be removed by a well planned program of surface dusting but most of it will require direct excavation. The yearly ablation loss is presently about 1 m a−1 of ice but Reference SlaughterSlaughter (1969) stated that as much as a 100% increase can be achieved by albedo modification.
Fig. 12. The calculated values of strain-rate, velocity and surface flux (or ablation) along line No. 6.
Discussion
The purpose of this study is to provide an estimate of the expected ice movement toward a proposed mining development. The configuration of the open pit is established at its final stage of operation and the flow of the glacier around the pit is calculated in order to provide an estimate of the volume of ice which must be excavated in order to keep the mine open. This study of the worst case of flow provides the necessary information for the feasibility study without calculating the flow at each stage of mine development.
The flow toward the pit is analysed by considering many segments around the pit and approximating the flow in each of these areas as two-dimensional. The equations are solved by averaging the stresses over the thickness of the glacier thus reducing the problem to the analysis of the average values of longitudinal stress, strain-rate and velocity. The basal shear stress is assumed to be constant over the entire glacier and certain longitudinal variations in the flow are neglected thus simplifying the task of relating the flow field to the configuration of the glacier.
The most serious weakness in this analysis is the close proximity of the reference position L, at which values of longitudinal stress and velocity must be assumed, to the margin of the glacier. The analysis given here is limited by the fact that topographic, ice thickness and velocity data have only been collected within a short distance of the deposit. Ideally the flow calculation should begin far enough away from the margin that the effect of the excavation would be negligible at that point. Further information of this type will be collected in the near future and then a better analysis of the flow at Isua can be made. The uncertainty about the basal temperatures in this area is also a matter of some concern since these have a direct influence on the parameters in the flow law; basal sliding (or even surging) is possible, and water discharge from beneath the glacier could flood the pit. Gudmandsen and Reference Gudmandsen and ChristensenChristansen (1972) reported that high temperatures of the ice at Isua prevented them from measuring ice thickness with radar, and the temperature measurements in bore holes suggest that the melting temperature exists at a depth of about 200 m. The ice along the margin is frozen to the bedrock. However, a short distance from the margin the basal temperature is probably at freezing point. The proposed margin for the full development of the mine will reach the 200 m depth between lines No. 4 and No. 5. Hence, further measurements of ice temperatures are required in order to locate the high temperatures more precisely. It is suggested that the area of ore labeled A in Figure 4 should be left in place in order to minimize the flow rates and crevassing near the margins. Clearly the problem with subglacial water provides a further incentive for sacrificing that ore. The development of the mine will take 30 years or more, a period during which much freezing of the basal water could take place. The presence of a large volume of water seems unlikely; the water probably exists in the form of a sheet which would take relatively little time to freeze. The temperature conditions at the base of the glacier will have to be monitored carefully during the development, since several competing factors will influence the basal temperatures: the snow cover will be removed more rapidly each summer, the ice thickness will be less, and the increased flow will increase the importance of both the advection of cold ice and generation of strain energy due to flow.
The ore body and associated rock formations are relatively resistant to erosion and thus provide a barrier to the flow of the ice sheet in this area. As the barrier is removed the flow will significantly increase due to changes in the surface slope and lowering of the margin. About 106 × 106 m3 of ice will have to be removed in order to expose the ore and an additional 66 × 106 m3 of ice will have to be removed in order to achieve the suggested profile. This profile will not be stable, however, and about 7.9 × 106 m3 of ice will have to be removed each year in order to prevent the glacier from thickening and advancing into the pit. Much of this ice can be removed by increasing ablation through surface dusting but continuous excavation will be necessary.
Acknowledgements
The Marcona Corporation through Mr Philip Morey, Supervisor of Mineral Resources Planning, sponsored this study. Mr B. Lyle Hansen and Mr Steven J. Mock made these arrangements possible and significantly contributed to the progress of the work. The computer program for the analysis was written by Mr Mock. Dr Malcolm Mellor reviewed this manuscript making many useful suggestions.
