The calving of icebergs from either tidewater glaciers or from floating ice shelves is the predominant ablation mechanism for existing ice sheets (Reference ReehReeh 1968, Reference Robin and deRobin 1979). It has been postulated that the stability of the marine West Antarctic ice sheet depends on the extent and stability of buttressing Ice shelves (Reference Ruddiman and MclntyreStuiver and others 1981). Disintegration histories of major ice sheets in the northern hemisphere indicate downdraw through ice streams and calving embayments as the significant mechanism of Ice-sheet removal (Reference Denton, Hughes, Denton and HughesDenton and Hughes 1981,Reference Ruddiman and MclntyreRuddiman and Mclntyre 1981). Ice shelves, because of their occurrence at high latitudes and general low relief, are considered more susceptible to slight climatic changes (Reference MercerMercer 1978). In this light, a better understanding of processes leading to calving is seen to be necessary for assessing both past histories of the cryosphere, as well as the future behavior of existing ice sheets.

Previous work has dealt with either analytic solutions for the configuration of isolated crevasses in an ice mass subjected to uniform tension (Reference WeertmanWeertman 1973) or an idealized ice shelf as a beam subjected to externally determined torques. Reference ReehReeh (1968) has studied the problem of a viscous beam subjected to a torque produced by the imbalance of hydrostatic pressures of the ice front, while Reference HoldsworthHoldsworth (1973)has modelled the stress and bending-moment distribution for a relatively Short ice ramp not in hydrostatic equilibrium due to a rapid change in water level. Reference IkenIken (1977) has used finite element analysis to describe a tidewater calving Situation.

In the present work, the calving process 1s investigated using a finite element program (Reference ZienkiewiczZienkiewicz 1971, Reference SchmidtSchmidt 1977) to obtain an approximate solution to the equations governing slow creeping Newtonian flow with velocity and pressure as independent variables. The program is based on composite quadrilateral elements formed by four triangular elements with the centroid of the quadrilateral a common node for the four triangles and the sides of the quadrilateral formed by the bases of the triangular elements. A linear velocity profile is assumed in each triangle with the hydrostatic pressure constant in the quadrilateral. The finite element technique used in the analysis presented here allows the specification of realistic ice configurations including cracks and spatially variable material properties. In addition one can use realistic equations of State (flow laws) without the constraints imposed by an analytic solution. Perhaps most important, one can look at the time evolution of the system as it responds to the changing stress configuration created by the development and propagation of a crack.

Of importance to note here is the difficulty that any of these methods encounter in attempting to quantify a fracture criterion. Stress and strain-rate distribution, as well as the past history of the ice in question, seem to affect the formation of crevasses (Reference HoldsworthHoldsworth 1969, Reference Hambrey and MüllerHambrey and Müller 1978) and estimates of the critical values for fracture range over Orders of magnitude. For this reason we have chosen the simplest criterion, i.e. that the crevasses form at the point of maximum tension. Propagation of the crevasse into the ice shelf will be discussed later.

Shown in Figure 1 is the finite element grid used to analyze a floating Ice shelf. The hatched area represents a highly viscous material with a density approximating that of ice-shelf ice, while the unhatched area represents the supporting sea-water in which the ice shelf floats. The densities of Ice-shelf ice and sea-water are taken to be 0.917 and 1.03 Mg m3 respectively. The flow-law constant of the ice-shelf ice is taken to be 3.14 bar a⁻¹ with the exponent equal to one, corresponding to an ice shelf 600 m thick at a temperature of -4°C (Reference ReehReeh 1968). The sea-water is represented by a Newtonian fluid with a flow law constant two Orders of magnitude lower than that used for the ice. While not the actual viscosity of water, this value produces accurate results while

Fig.1.

Finite element grid used to model a floating ice shelf. The unhatched region is a low viscosity water-like material while the hatched region is a higher viscosity ice-like material. The point marked G.L. represents the grounding line and is constrained not to move.

avoiding a round-off error problem due to the actual difference in order of magnitude of 15 in viscosities of ice and water. The lower left corner of the ice-like material (marked G.L.) represents the grounding line and is constrained not to move (a finite grounding-line velocity in the right-hand direction would not affect the overall analysis). An important aspect of this model is that because both the ice and water are included, the boundary conditions at the icewater interface are automatically included in the analysis. The initial configuration is one in which hydrostatic equilibrium exists at every point along the flow line. Boundary conditions for the model then consist simply of a free top surface, no flow into or out of the water region (right side and bottom) and the left-side boundary constrained only to move in the vertical direction (no flow across the grounding line). Due to the horizontal imbalance of hydrostatic forces at the ice front a bending torque is generated. Unlike Reeh's model in which this torque is externally calculated and applied as a boundary condition, this bending torque is generated internal to the model by the interaction between the ice elements and the water elements in the finite element program. Figure 2 shows a comparison of deflections from Reeh's beam analysis with deflections of the centerline (midpoint from top surface to bottom surface to compensate for ice thinning allowed in this model but not dealt with by Reeh) predicted by the finite element method. Qualitatively they predict very similar results. The bending down of the ice front is compensated for by the arching up of the ice surface from one to several ice thicknesses back from the ice front. This arching and bending produces a non-uniform distribution of tensions along the top surface with a maximum tension of 3.44 bars about one half an ice thickness or 300 m from the ice front. This tension at the surface is shown in Figure 3. Continued time sequencing (using the output velocities to distort the grid which is then used as input for the next time step) shows little change in this stress distribution. Assuming that a crack will form at the Position of maximum tensile stress, we introduce a crack by removing an element 10 m wide at the point of maximum tensile stress. Shown in Figure 3 are curves of surface tensile stress for cracks 30 and 60 m deep. Tensile stress is reduced on either side of the crack while tensile Stresses more than one ice thickness away are not appreciably affected. This leaves the possibility of a second crack developing about one half an ice thickness behind the first crack at the new point of maximum tensile stress. We will not deal further with this possibility. Time evolution of the surface tensile stress after the crack has opened is shown in Figure 4 for a crack with bottom at sea-level (65 m deep). Little change is seen in the vicinity of the crack, the only notable effect being a slight increase in the tensile stress far from the ice front.

Fig.2.

Comparison of deflections after 0.325 a has elapsed predicted by the viscous beam approximation of Reeh and by the finite element method.

Fig.3.

Normal stress at the centroids of elements along the surface of floating ice for no-crack and for cracks of different depths. The grounding line is to the left and the calving front to the right.

Fig.4.

Normal stress in the surface elements from time of 0.01 to 0.03 a for a crack extending to sea-level.

Lacking an adequate fracture criterion for ice, we were unable to automate the propagation of the crack into the ice. Instead a heuristic approach was taken whereby two independent criteria were applied to determine if the crack would increase in depth. The first criterion was that the crack be opening along its entire depth. This was necessary because it was found that a crack of up to 115 m deep would still be opening at the bottom while at an intermediate level (slightly above sea-level) it would be closing. The opening at greater depths is thought to be caused by the general thinning and elongation of the ice shelf, an effect not dealt with by Reference WeertmanWeertman (1973). Weertinan's prediction for the depth of a water-free crack in an ice mass subjeeted to a 3 bar uniform tension (the average of the tension produced by the unbalanced hydrostatic forces at the ice front of an ice shelf 600 m thick) was 52 m, which agrees well with our prediction of a stable, non-closing crack with bottom close to sea-level (about 60 m deep).

The second criterion adopted to determine continued crack propagation was that tension be observed in the ice element immediately below the crack bottom. The finite element prograin Outputs reliable stress and strain-rate Information at the centroid and at the four integration points within each element. Thus for an element 30 m deep below the crack one can reliably determine the tensions at three levels within the element, 6, 15, and 24 m below the crack bottom.

As has been pointed out, a water-free crack can only open to approximately sea level in an ice shelf 600 m thick. Reference WeertmanWeertman (1973) provides an analytic expression for the depth of propagation of a crack filled to various depths with water (the hydrostatic pressure provided by the higher density water serves to balance the hydrostatic pressure in the ice). For a crack filled two-thirds full of water, the analytic solution yields a crack depth of 126 m. Applying our two criteria to a crack two-thirds full of water (done in the finite element program by applying calculated hydrostatic pressures to the faces of the elements surrounding the crack) one finds the crack stable and opening to a depth of 138 m. Applying the criterion of necessary tension in the element below the crack, one finds the level of zero tension coineiding with the crack bottom at a crack depth of 113 m. We accept these results as relative quantitative agreement with the analytic Solutions.

Dp to this point little has been done that could not be performed equally well by the analytic methods. The strength of the finite element method lies in continuing the analysis in time beyond the point where the assumptions necessary for the analytic Solutions have broken down. Shown in Figure 5 is the time evolution of the configuration of the crack with the left-hand edge of the crack repositioned to zero to remove the displacement down-stream produced by thinning of the ice shelf. The time step between adjacent configurations is approximately 2 d. In Figure 6 one sees the corresponding tensions at different levels within the element directly below the crack bottom. One can readily note a general trend of increasing tension in the element as time progresses. If one assumes a linear Variation of stresses within the element one can obtain a measure of the depth of the level of zero tension. The movement of this level of zero tension with time is shown in Figure 7 as the relative level change, plotted three ways, with respect to sea-level, with respect to the crack bottom, and with respect to the ice surface. The various curves have been displaced for clarity. From these curves rates of movement of the level of zero tension can be obtained. For this particular crack depth of 108 m one finds the level of zero tension moving deeper into the ice at a rate of between 0.41 and 0.30 m d⁻¹ depending on the reference point chosen. A similar analysis for cracks of different depths shows similar lowerings of the level of zero tension with time, shown in Figure 8. This analysis suggests that, while absolute tensions below an opening crack may decrease to zero as the crack propagates deeper into the ice, the time evolution of the ice shelf will tend to increase the tension below a crack. Thus, a crack which extends to the level of zero tension at a certain time, will, as time passes, develop greater and greater tensions immediately below the crack. In the absence of a well-defined fracture criterion one can hypothesize that the crack might propagate rapidly to near the level of zero tension and stop until the ice shelf as a whole has deformed with time to a sufficient extent to again initiate fracture and crack propagation. This deformation leading to increased tension below the crack appears to accelerate as the crack becomes progressively deeper .

Fig.5.

Time evolution of a crack, initially 10 m wide and 44 m deep, filled two-thirds full of water, from time of 0.015 to 0.05 a.

Fig.6.

Normal stress with time at various levels within the element directly beneath the crack bottom. The straight lines represent least Squares fits to the data.

Fig.7.

Movement of the level of zero tension as a function of time taken with respect to sea-level, the crack bottom, or the surface. Least Squares fits are shown by straight lines

Fig.8.

Rate of movement of the level of zero tension for cracks two-thirds full of water and differing bottom depth below sea-level. Zero-level movement is measured with respect to surface, sea-level, and crack bottom.

As an extreme example of this phenomenon we took a crack 530 m deep (about 90% of the ice thickness) filled with sufficient water (90% of the crack depth)so that the crack bottom was near the level of zero tension. Time evolution of this configuration showed a very rapid increase of the tensions in the element directly below the crack with the level of zero tension moving as rapidly as 10 m d⁻¹.

With this Information one can hypothesize a reasonable scenario for crack evolution leading to calving. Following the calving of an iceberg the ice front will deflect downward at the same time that the ice surface more than one ice thickness from the ice front arches upward. This will lead to a tension maximum about one-half to one ice thickness from the ice front. Fracture is assumed to occur in this region. In the absence of melt water this crack can only propagate to near sea-level. Indeed, throughout the entire evolution of the crack it will have to be filled to near sea-level to prevent the crack from closing. It has been seen that the time evolution of the distorting ice shelf will somewhat relax this condition as the ice shelf bends in such a way as to increase the tension at the crack bottom with time. Iceberg calving by this mechanism will be strongly seasonal in character since it depends on copious melt water to keep the crack open. Not dealt with in this analysis is the possibility of bottom crevasses, which can extend through close to 78% of the shelf thickness (Reference WeertmanWeertman 1980) joining with surface cracks to produce calving.