I. Introduction
Many floating ice masses show strikingly regular patterns of melt pools during the summer ablation season. Figure I shows the pools on George VI Ice Shelf, Antarctica (lat. 71°S., long. 68° W.) ; similar pools have been reported on the Ward Hunt Ice Shelf, northern Canada {Reference Koenig, Koenig, Greenaway, Dunbar and Hattersley-SmithKoenig and others, 1952, p. 66) and Bessels Fjord, north-east Greenland (Reference Post and La ChapellePost and LaChapelle, 1971, p. 101). The pools apparently form in the troughs of regular surface undulations, of wavelengths 450 m, 250 m, and 40 m for these three cases, and although many explanations have been advanced for the formation of these undulations (Hattersley-Smith, 1957), none can be regarded as conclusive. While searching for a more satisfactory solution, an analysis has been made of thermal stresses set up in the top 5 or 10 in of an ice shelf or indeed any glacier as a result of thermal expansion and contraction due to surface temperature variations. Temperature variations due both to the annual temperature cycle and to typical cyclonic conditions have been taken into account. Such temperature fluctuations have been mentioned before in connection with surface undulations on ice shelves (Reference ZubovZubov, 1955, p. 41) and on sea ice (Reference Wright and PriestleyWright and Priestley, 1922, p. 344), and thermal cracks in sea ice have been explained using the model of flexure of a thin elastic floating plate (Reference Evans and UntersteinerEvans and Untersteiner, 1971). Here we are concerned with ice of great thickness. Whilst it appears from the following analysis that thermal stresses cannot be invoked as an explanation of periodic surface rumples of the wavelength found, the treatment and results reveal that thermal stresses are dominant in the top 2 or 3 m of a glacier and may result either in the formation of shorter-wavelength rumples or in fracture under tension.
Fig. 1. Melt pools on George VI Ice Shelf, Antarctica, at separation approximately 450 m. (Photograph by Charles Swithinbank.)
2. Temperature Variations Near the Surface of a Glacier
We shall first look at the way in which temperature variations at the surface of a glacier are transmitted downwards by heat conduction. This is a standard problem in heat conduction theory: we begin by taking a simple sinusoidal variation of surface temperature and apply this as a boundary condition to find the appropriate solution to the Fourier equation of heat conduction (Reference Carslaw and JaegerCarslaw and Jaeger, 1959, p. 64-66; Paterson, 1969, p. 169). Let us begin by taking a surface temperature of:
where T is temperature as a function of depths and time i, T0 is the mean annual temperature, TV is the amplitude of the temperature variation and cur/27r is the frequency of the variation.
We require solutions to the Fourier equation:
where k is the thermal diffusivily of ice.
Neglecting transient time-dependent terms in the general solution, which are important only when the temperature fluctuation is first set up, we find the steady-state solution of Equation (2) for boundary conditions given by Equation (1) to be:
Solutions of this form for different phase, frequency, and amplitude are superposable, so we can construct a solution for any surface boundary condition which can be expressed as a series :
where is a dimensionless phase term. This is equivalent to saying that any boundary condition can be treated which satisfies Dirichlet's conditions and is therefore susceptible to Fourier analysis; this will of course always be true for any realistic surface temperature conditions provided we allow sufficiently many terms in the expansion, but the handling is much simplified if we can idealize the real boundary condition into a sequence of only a small number of terms in Equation (4).
Two things should be noticed about the character of the solution given by Equation (3). First, any surface variation of temperature becomes attenuated with depth according to the factor
The presence of a — (cor)k term in the exponent means that rapid temperature fluctuations are more quickly attenuated than slower ones. Because of this we find that diurnal temperature variations are attenuated to 5% of their surface amplitude by the time they reach 50 cm depth, while the annual cycle is attenuated to 5% only when it reaches zo m depth. It is therefore found that at most depths the annual cycle dominates other shorter-period variations. Secondly, the phase term(W712k)i has the effect that more rapid fluctuations are more rapidly transmitted; thus we find that the temperature minima of a variation with a period of sixteen days propagate to 3 m depth in a time of 11.3 d, while the minima of the annual cycle take 53 d.
3. Temperature Data and a Simple Model
We shall first carry out an analysis for the conditions prevailing in the area of George VI Ice Shelf where the regular patterns of melt pools are observed, and later consider alternative situations. Actual surface temperature data may be found from synoptic meteorological observation results for air temperature, which are available for the Fossil Bluff region of George VI Ice Shelf throughout the year 1968 (personal communication from D. W. S. Limbert). The mean daily temperatures are plotted for the year in Figure 2(a). Superimposed on this there is also the continuing diurnal variation, which is of amplitude about 5 deg, but this has been omitted since such high-frequency temperature variations are effectively attenuated within the top metre of ice. Any calculations concerning ice very near the surface are in any case unjustified because of the unknown effects of air percolation, radiation absorption, melt percolation, and the drifting of surface snow. Also displayed (Fig. 3) is a graph of half-monthly mean temperatures, with an annual sinusoidal curve fitted to them. The fitted curve has an amplitude of 11 deg about an annual mean of IOC, or 262 K. The plot of daily temperatures throughout the year shows that this simple annual sinusoidal model is a gross simplification, for in fact quite long and marked cold and warm spells occur throughout the year, especially in winter, and there is a large day-to-day variation within these spells. The amplitude of these shorter-period variations is of the same order as that of the annual cycle so that they cannot a priori be neglected.
Fig. 2. Mean daily temperatures throughout the year: (a) data recorded at Fossil Bluff, 1968; (b) empirical model reproducing principal features {Equation (5)). The dashed line shows the mean annual temperature.
Fig. 3. Half-monthly mean temperatures for Fossil Bluff, rg68 displayed with a fitted annual sinusoidal curve and a dashed line to show mean annual temperature.
An empirical model has been constructed to take better account of the most important features of the data. Such features are:
(i) the annual cycle, of amplitude it deg;
(ii) an irregular variation of period in the region of two months, amplitude some 7 deg.This shows up in Figure 3 ;
(iii) an irregular cycle of period about 17 d, amplitude about 7 deg; (iv) smaller variations over periods of between 3 and 7 d, of amplitude about 3 deg.
We shall not enter here into discussion of the origin of these quasi-regular weather cycles, which are presumably associated with the natural periods of stability of cyclones in the region, since it is enough for our purposes to know that they exist. That they are real features can be seen from a spectral analysis of the data. Such an analysis has been carried out using the technique of maximum entropy, and the power spectrum (Fig. 4), shown here considerably-smoothed to reduce noise, displays significant peaks at periods of 16.7 d, 6.2 d, 3.9 d, and 3 d. The variation (ii) of period about two months, could not be resolved with a single year's data, but is included in order to simulate the marked warm spell, the "kernlose winter", which is a characteristic feature of Antarctic weather patterns (Van Loon, 1967). It is apparent here between June and September in Figures 2(a) and 3.
Fig. 4. Maximum entropy power spectrum of mean daily temperature variations for Fossil Bluff, 1968, showing power plotted logarithmically against period.
A further important feature is that all the subsidiary fluctuations (ii), (iii), and (iv) are more pronounced in winter than in summer, which is due to the resistance afforded by the large mass of ice in a polar environment to any tendency for the surface air temperature to rise much above o°C ; it acts as a buffer. This characteristic has been included by means of a modulating envelope of half period one year, which multiplies sinusoidal variations of period 62 d, 16.7 d, 6.2 d, and 3.9 d. The full expression is then:
with r0= -II°G = 262 K,
The amplitudes of the various waves have been chosen simply by eye and trial, and the resulting curve is shown in Figure 2(b). The agreement in general form, although not in particular detail, is readily seen. It is sufficient that the model represent a possible year's temperature data; there would be no great surprise if one year they occurred according to the model given.
We are now in a position to find the form of temperature with depth for Equation (5) applied as a boundary condition. To do this we need to force Equation (5) as nearly as possible into the form of a series as in Equation (4), which is then readily treatable. This can be achieved by noticing that 1/2ω, is at least an order of magnitude smaller than any of ω2, ω}, αι4, or ω5, and that we can therefore treat the envelope sin (w1t/2) as a slowly varying function of t. This means that we treat it as a constant for the purposes of differentiation with respect to t and so can write Equation (5) as:
where we regard terms Ts sin (ω,//2) as the constant terms Tr in the expansion in Equation (4).
It is then straightforward to proceed directly to the complete steady-state solution as a function of depth and time. It is:
where k is equal to 1.091 X 10-6 m2 s_I. This is plotted in Figure 5 to show temperature variation over the year at depths of 0.5, 1, 2, 3, 4, 5, and 10 m. The effects of selective attenuation and phase delay can be seen clearly.
The value used in Equation (6) for the thermal diffusivity k is calculated for ice typical of George VI Ice Shelf. Diffusivity is strongly density dependent and this dependence can be expressed by the empirical formula (Van Dusen, 1929):
where c is the specific heat capacity of ice, taken as 2.01 × 103 J kg—1, K—1, and ρ is the density, expressed in Mg m—3. Density profiles have been taken for George VI Ice Shelf to 10 m depth, and show the ice in the melt area to be of quite uniform high density, lying within the range 0.90±0 0· 2 Mg m—3 (personal communication from J. L. W. Walton). This very high value for surface density is a result of the extensive summer melt which characterizes most of the region: almost all the year's snowfall, some 50 cm, is melted each summer and either runs off the ice shelf down tide cracks and moulins or forms pools on the surface. Any net accumulation over the year, which seldom amounts to more than a few centimetres (personal communication from J. F. Bishop), is a result of refreezing of melted snow. It is interesting to note that these conditions correspond closely to those on the Ward Hunt Ice Shelf (Hattersley-Smith and others, 1955, p. 22 23). Under the more commonly found polar conditions of snow and firn in the top ten metres we expect lower values for the purely conductive thermal diffusivity (k = 0.613 × 10–6 m2 s1 for ρ = 0.6 Mg m-3; k =0.407 × 10–6 m2 s-1 for ρ = 0.4 Mg m-3) ; however, it is found from field measurements on the rate of propagation of cold waves through snow and firn (MacDowall, 1964) that the presence of air passages enables convective processes of propagation to compensate for the reduced conductivity of low-density snow. The effective diffusivity is almost the same as that for pure ice, if not slightly greater. Our results for (temperature-wave propagation should therefore be applicable to all types of surface cover.
Fig. 5. Progress of a model temperature wave through the top 10 m of a glacier.
For the purposes of calculating thermal stresses due to rates of expansion and contraction of the ice mass as its temperature rises and falls we shall be concerned primarily with the temperature-time gradient dT/dt as a function of time and depth. This can be calculated directly from Equation (6) by differentiation, and we shall again neglect the small terms arising from the differentiation of the envelope term sin (ω,t/2). We find:
An important feature of the form [inline img 1] with depth is that where terms in the expression for T(y, t) attenuate according to exp [—y(ωs/2k)½] terms in [inline img 2] attenuate according to o>sexp [—y(ωs/2k)½] which is a product of two terms in ωs, one increasing and the other decreasing. So, although the exponential term will always eventually dominate the linear term, the temperature-time gradients due to quite short-period fluctuations can still be important at a or 3 m depth; it is not the case, as it is with T(y, t), that the longest-period wave dominates at all depths. We can in fact find the frequency of surface temperature variation which for a given amplitude gives the maximum values of [inline img 3] at a given depth: by maximizing the amplitude term T sωs exp [—y (ωs/2k)½] with respect to ωs we get (ωs = 8k/y 2 for the value of ωs giving maximum [inline img 3] at depth y. At 1 m depth the wave of period 8.5 d dominates; at 2 m, 33 d; at 3 in, 75 d; and at 6 m, good.
The treatment so far rests heavily on the assumption that a conductive model is valid; it has taken no account of heat transfer due to melt-water percolation and refreezing. This has been found to be extremely important on, for instance, Isachsen's plateau, Spitsbergen (Sverdrup, 1935), where the percolation of summer melt through firn caused a very rapid elimination of the winter cold wave, amounting to values for [inline img 4] at 3 m of I deg d_I over a period of about a day and of 0.5 deg d–1 over a period of about ten days. These compare with values using Equation (7) of 0.1 deg d -1 and 0.07 deg d–1 at 3 m for similar periods. However, for George VI Ice Shelf it is clear that we are not in a region of permeable firn, so that the contribution of percolating melt water is likely to be very small; further, since any contribution there is will tend to accelerate the warming process, we can regard the calculations of stress which follow as setting a lower limit on those actually present.
4. Thermal Stresses in High-Density Ice
Theory of thermal stresses
The problem we wish to solve is that of the stresses in a thermally fluctuating region of some 5 m thickness attached to a large mass of ice which is thermally stable but which may anyway be creeping longitudinally under the stresses of external constraints and its own weight.
Let us begin by looking at the simpler case of a thermally active region at the surface of a completely static, unstraining larger block of ice. Because the surface layer is much greater in lateral extent than in thickness and is everywhere joined to the larger mass (indeed they have no well-defined boundary), we can say that any temperature changes in the surface layer cannot give rise to large-scale movement (expansion or contraction) of the layer relative to the stable block beneath. This is an application of Saint Venant's principle concerning edge effects: if a surface rests on and adheres to a massive block as described and undergoes thermal variations, we should expect the tendency towards thermal expansion and contraction to lead to actual motion only in the vicinity of a free boundary. Far away from the edges no relative motion of the surface layer and the block can occur. The same condition of zero actual strain applies in the case of sea ice or lake ice which is restrained by rigid banks at its boundary. We are therefore treating the case of laterally restrained thermal expansion and contraction: any tendency of the ice to expand or contract results not in strain but in the creation of just enough stress to restrain the tendency. It will also be assumed that the ice does not fracture under tension, an assumption which will be justified once we have calculated the stresses involved.
In order to calculate the stresses we use a method in which we in effect transform our spatial frame of reference to an expanding frame of coordinates. We choose these expanding coordinates so that they expand at the same rate as would an unrestricted body undergoing thermal expansion: they therefore represent the frame in which the material in question would be stress-free. So, if the material is subject to a temperature-time gradient of [inline img 5] and has a linear coefficient of expansion a, it tends to remain stress free by expanding at rate [inline img 6] To treat it we therefore transform from real-space coordinates Xi to expanding coordinates Xi related by:
If we then require that the material remain unstrained in real space we are introducing an effective strain-rale of:
relative to the stress-free frame xi and are therefore introducing the stress associated with it.
We shall assume that this stress is simply that required for the strain-rate concerned according to the empirically determined flow law for ice. This assumption is not exact, although it is one generally made in related literature (for instance Gatewood, 1957, p. 1). An exact analysis involves taking complete account of the thermodynamics of irreversible processes, and finding values for Onsager's "cross-coupling" coefficients for a material with a non-Newtonian flow law (Boley and Weiner, 1960, p. 1-40). If such an analysis were possible it should lead only to minor modifications of the stresses involved (Boley and Weiner, 1960, p. 41-44); the "cross-coupling" is an expression of a relationship analogous to the complementary relationship between stretching and heating a rubber band: if much expanded it will warm slightly, and if much warmed it will contract slightly. Only very tiny errors are involved, however, in calculating the stresses in a stretching rubber band on the assumption that it remains isothermal.
The advantage of the formulation in terms of the expanding coordinates is that we can now readily treat the real case in which the thermally active layer rests on top of a continuum which is itself undergoing gross strain. If the ice mass as a whole is undergoing longitudinal strain-rates ⋵i relative to the real coordinates xi, then it is undergoing longitudinal strain-rates
relative to the thermal stress-free system x/. We shall not enter here into any treatment of shear stresses: the shears in the vertical plane are in the case of an ice shelf anyway negligible, and shears in the horizontal plane may be eliminated from the analysis by choosing xi along the principal axes. The orientation of the principal axes is unaltered by the transition to the frame xi
We are now in a position to make an analysis of stress. An elastic model is unrealistic, since ice responds elastically only over a time scale of 5 to 10 s (Gold, 1958). We shall therefore consider creep, and for high-density ice we shall use as flow law the form most widely accepted on the basis of theoretical expectation and experimental confirmation, that is (Paterson, 1969):
In this equation, expressed here relative to the system AT/, we have an expression for the strain-rate $$$$i as a function of the stress deviator $$ oj*, the temperature T, and the effective stress T. The quantities A, Q, and n are empirically determined constants and R is the gas constant.
The use of the stress deviator $$at* instead of the stress σί expresses the observed fact that strain-rate is independent of hydrostatic pressure, and is related only to the deviation of stress from hydrostatic pressure (Rigsby, 1958). The two quantities are related by:
The Boltzmann exponential term in Equation (9) expresses the increasing ductility of ice as temperature increases (Glen, 1955).
The effective stress τ arises from the second invariant of the stress-deviator tensor, which is expressed by:
The term τη * appears in the flow law as a result of general requirements on the form of flow law allowable (Nye, 1957; Glen, 1958): because randomly oriented polycrystalline ice (which is what we are assuming glacier ice to be) is isotropic the law must be independent of axis orientation and therefore be expressible as a general relationship involving only the three invariants of the stress deviator and the strain-rate tensors. On requiring that under given stress conditions the components of strain-rate be proportional to the components of stress deviator, following the Lévy-Mises equation, we find that we may eliminate terms involving the third invariant; and we find also that by Equation (10) the first invariant vanishes identically, so that we are left with a relationship involving only the second invariants. We posit :
with ê' arising from the second invariant of the strain-rate tensor
and it has been found by experiment that n lies in the range between 1.5 and 4.2 for conditions of steady-state deformation. The experimental results are summarized in, for instance, Weertman (1973). We shall ignore here the possibility of fatigue hardening, for which no evidence exists, and shall ignore any transient flow characteristics ; such transients correspond to the time taken to establish stable deformation mechanisms in secondary creep and generally die out in a few hours (Glen, 1955). We shall use here the value of n = 3.
For the purposes of the analysis in hand we can treat stresses as invariant under the transformation between frames xi and xi. This is justified because stress has dimensions
Now mass and time are invariant between (slowly moving) inertial frames, and although length is not strictly invariant since we are transforming between relatively expanding frames, it will be seen that we are dealing with strains in the order of 10-3, so that errors arising from the treatment of stress as invariant are tiny.
Because the volume of ice is changing we no longer have the familar equation for the null divergence of the velocity field
but instead have the equation
where 3a dT/dt represents the rate of change of volume per unit volume. However, by substituting Equation (8) into Equation (13) we find
in the stress-free frame, as one might expect.
Through Equation (9) this means we have also:
We proceed as follows. At a point, the factor A exp [- Q/RT] τζ$$$ in Equation (9) is the same for all components $$ aj*, and we can therefore write it as y and say :
and
Hence
or
where
From Equation (14) we can write
and hence by Equation (15)
so that we can express τ in terms of a and σ,*. Shear stresses are being assumed negligible so that Equation (11) becomes
and therefore we have
by virtue of Equations (15) and (16), leading to
Our expression for strain-rate as a function of stress deviator is then, by Equations (9) and(17):
We can therefore calculate thermal stresses according to the formula
and use Equation (15) for the other principal stress a2*
We notice that ι+α+α2 is necessarily positive non-zero, and so cannot lead to a singularity.
Numerical results for stresses
Equations (18) and (15) can now be used to calculate the thermal stresses associated with realistic values for the parameters T, dT/dt, iu and e2. To begin with we shall look at the case of an ice sheet not itself undergoing gross strain. That is :
and therefore by Equation (8)
Thus
and the stresses are
with a = 5.3 X 10-5 deg-' (Hobbs, 1974, p. 347), Q/R = 7.293 X 101 deg (Weertman, 1973), and A = 1.05 x 1o-12 N-3 m6 s-1 (from Paterson, 1969, p. 83.
Fig. 6. Progress of the stress wave-form though the top 10 m ofa glacier, fir the case of a glacier not undergoing gross strain.
The value used here for A is at the upper end of the range found experimentally (Weertman, 1973); stresses calculated from it are therefore the minimum to be expected. Figure 6 shows the stress as a function of time of year for depths 0.5, 1, 2, 3, 4, 5, and 10 m. Figure 7 shows the values of (a) the maximum magnitude of stress $$ lc*|max occurring over the year, (b) the mean magnitude of stress $$|σ*| over the year and (c) the overburden pressure — pgy due to overlying ice. It is seen that in the top 3 m of surface cover, stresses of thermal origin dominate those due to the weight of ice above. The mean stress over the year taking account of sign is very nearly zero: it is not precisely zero, since the symmetry of the surface temperature form is destroyed by the effects of selective phase delay of the temperature waves, and this combined with the cubic form of the flow law leads to a net bias of compressive stress over extensive stress at most depths. The magnitude of the effect is however entirely negligible, having a value at 3 m of 0.005 bars.
Fig. 7. Comparison of thermal stress with overburden pressure, as a function of depth.
We consider next the case of an ice shelf undergoing gross strain. Strain-rates measured on George VI Ice Shelf are typically in the order
(personal communication from J. L. W. Walton). The minus sign denotes compression.
Use of these parameters in Equations (18) and (15) leads to the graphs in Figure 8. These results show that in the top 2 or 3 m of a glacier the thermal effects lead to stresses alternating quite rapidly between extensive and compressive, and these alternations are enough to dominate the stress due to dynamic creep even when there is an overall strain of some ±10 10 s r present. Such a strain is moreover at the upper end of the range of strain-rates to be expected on ice shelves (see for instance Thomas, 1973, Table I).
We can therefore safely say that, in the top 3 m of ice cover, stresses due to temperature change dominate those originating both in gross deformation and in weight of overlying ice.
Fig. 8. Progress of the stress waveform through the top 10 m of a glacier, for the case of a glacier undergoing gross strain-rates $$ il — ~1.9 ° s~' and έζ =0.6X io~10 s~:. The solid line shows the stress along the major principal axis, and the dashed line shows the stress along the minor axis
Effect of the stresses
We now have a picture of the surface cover of a glacier undergoing quite pronounced alternating compressive and extensive stress throughout the year. What effects can such stresses have? First of all let us determine whether our assumption that no tensile fracture occurs was justified, and then proceed to consider the possibility of other modes of permanent deformation.
We have found that thermal stresses in high-density ice are generally in the order of 0.5-1.0 bar, and in any case always less than 4 bar; the corresponding strain-rates are less than 5X 10-9 s_I. The fracture behaviour of ice has been studied under conditions both of constant applied stress and of constant applied strain-rates. Butkovich (1954) finds the critical fracture stress under tension to lie in the range 13.7 to 18.7 bar for the temperature range 0 to — 40°C, and Gold ([1971]) finds that the critical strain-rale under compression is about 10-7S-1. We can therefore say that for high-density ice in the temperature range considered it is unlikely that fracture under thermal stress can occur unless local inhomo-geneities cause stress concentrations.
The assumption of no fracture is therefore validated and we can consider the effect of the stresses as calculated. Let us look at the possibility of the formation of surface rumples due to buckling and boudinage under stress. This would be an effect analogous to the formation of rumples in the surface tarmac of a road during a heat wave.
It might at first be objected that the magnitude of strain involved, which is in the range 10-4 to 10-3 in the surface four metres, is too small to result in fold formation of significant amplitude, for geological orogenesis is generally concerned with continuing compressive strains of the order of 5% to 25%, or about 10-3. A strain of 10-3, which anyway fluctuates about zero strain, might seem unlikely to result in folding of appreciable magnitude. An order-of-magnitude calculation shows that in fact such strains need not be insignificant: consider the case of a flexible incompressible bar, pinned at its ends and undergoing thermal expansion (this is the well-known problem of the buckling of railway tracks when warmed). Assuming that a straight bar distorts into an arc of a circle we find (Fig. 9) that a small strain e in a bar of length L results in a deflection h at the highest point of
Fig. 9. Buckling of a pinned incompressible bar subject to compressive strain.
For the case of strain 10-3 this amounts to o.o19L, or a deflection of 8.7 in for rumples of wavelength 450 m, and of 19 cm for rumples of wavelength 10 m. These are large deflections. The model of a pinned incompressible bar does not, of course, correspond closely to reality; nor does the model we used above to calculate stresses on the assumption that deformation is perfectly homogeneous pure shear, for it emerges that such deformation is generally unstable to the formation of folds. Reality lies somewhere between the two models, so that we should expect some buckling inhomogencity to occur, but with an amplitude less than 0.019L; the model of a flexible bar sets an upper limit, perhaps a generous upper limit, on the allowable amplitude.
Let us try to calculate a wavelength for buckle formation. There is a fairly extensive literature on buckling and orogenesis of visco-elastic media under stress, due largely to Biot (1961, 1965). The resulting predictions have been used before for glaciological applications: the surface undulations on the Ross Ice Shelf have been modelled as a gross folding of the ice shelf under compression (Kehle, 1964), and surface rumples on the Meserve Glacier, Antarctica, have been modelled as due to compression of a medium with exponentially decaying viscosity with depth (Holdsworth, unpublished). More recent work (Smith, 1975, 1977) has stressed the importance of the fact that the material in question is non-Newtonian and it is Smith's development to which we shall refer.
We shall need to look at the profiles of stress and viscosity with depth at various times during the year. For the "viscosity" of a non-Newtonian material we shall adopt the parameter "effective viscosity", which characterizes the relationship between an increment in stress, $$ δσ*, and the resulting increment δέ in strain rate. For an increment in longitudinal stress acting against a background of longitudinal stress this quantity is given by
or, using the results found earlier in this Section,
and, for a = 1,
Fig. 10. Stress as a function of depth for six days of the year.
Inspecting the profiles of stress and of viscosity against depth (Figs 10 and 11) we see that the form is essentially that of a series of two or three layers of some 2 m thickness, alternately in compression and extension at stresses between 0.5 and 1 bar. The form of Equation (19) used for the viscosity leads to singularities as the stress deviator σ,* passes through zero, implying that the viscosity becomes infinite at these points; this feature is due to the cubic form adopted for the power flow law, and disappears if we consider ice to show more nearly Newtonian behaviour at low stresses, as found by Butkovich and Landauer (1960). Apart from this feature, the viscosities of successive layers show a steady but quite minor increase with depth.
Fig. 11. Effective viscosity as a function of depth for six days of the year.
Figure 12 shows the type of stress model we would wish to analyse, although, since the form of stress with depth is continually changing, we cannot expect to make an analysis of buckling behaviour which will be precisely valid throughout the year. The model of Figure 12 is further complicated by the fact that our analysis of stress has considered successive layers to be in differently expanding coordinate frames, so that there are difficulties in expressing the boundary conditions at an interface. The treatment of buckling behaviour proceeds by considering the instability of small perturbations in the interface between two layers; we find out the rate at which a small perturbation grows or decays under the stress conditions applying, and thus determine the initial growth rate as a function of the wavelength of the perturbation. The wavelength which has the greatest initial growth rate is identified as the "dominant wavelength", and it is this wavelength which we should expect to find occurring. From the analysis of Smith (1975) we can see that the primary driving force of interfacial instability arises from a discontinuity in longitudinal stress deviator in passing from one layer to another; in his analysis, which is for homogeneous pure shear of a stratified viscous block, the stress discontinuity arises through the requirement that a stratified block of layers of different viscosities deform at a uniform strain-rate. The layers must therefore sustain different longitudinal stress deviators. In our case we again have an interfacial stress discontinuity, indeed a larger one, in which successive layers are under stresses of opposite sense. In the absence of a precise analysis, the most we can say is that such a physical situation would seem very likely to lead to buckling instability. As to the wavelength, we can make a fairly firm order-of-magnitude estimate by making a general survey of the results from the many possible models of buckling which have proved susceptible to an exact analysis (Fig. 13). It is found in general that the dominant wavelength is primarily a function of the layer thickness, and is very insensitive to other factors: it depends in a minor way on the viscosity ratios of successive layers, on the conditions of stick or slip at the boundaries, and, if gravitational forces are taken into account, on the ratio of deviatoric stress to gravitational stress. It is found, however, that for all but the most extreme cases we have a ratio L/d of buckling wavelength L to layer thickness d of between 4 and 6.
Fig. 12. Schematic diagram of a typical surface thermal stress state.
Since we are concerned here with layers of thickness approximately 2 m we can assert that any surface rumples which are the result of thermal stress may be expected to be of wavelength between 8 and 12 m. Thermal stresses do not, therefore, provide an explanation of the rumples of the wavelength found on ice shelves.
Analyses of buckling behaviour invariably restrict themselves to "incipient" buckling-— that is, the initial growth rate of infinitesimal perturbations. No attempt is made to quantify post-buckling behaviour, and for this reason no estimate is made of the amplitude of folding at a late stage in the process. The fact that in our case the ice is undergoing periodic extension and compression instead of steady stress may increase the growth rate of rumples: for the formation of rumples under compression will result in regions of thinner ice in which stress will concentrate when the system is placed under tensile stress. The rumples should be accentuated by a process of necking. There is also the possibility of melt water and accumulation finding their way into fissures in the ice (Shumskiy, 1955) ; such processes would render the process of expansion and contraction irreversible, and lead to compressive stress dominating. This is especially true in the case of melt-water percolation, which will lead to additional stresses due to the volume change on freezing.
For an unstraining ice mass we can say nothing about the preferred orientation of rumples since there is perfect symmetry in the horizontal plane, but we might expect quasi-hexagonal patterns as with cracks in frozen ground. For a straining ice mass we should however expect the rumples to form perpendicular to the principal axes of longitudinal stress, for these are the directions in which the total stresses are a maximum.
Fig. 13. Modes of buckling under stress: (a) viscous layer resting on a semi-infinite viscous medium (Biot,1961); (b)non-Newtonian layer embedded in an infinite non-Newtonian medium (Smith, 1975,1977). The shaded areas represent the more competent (more viscous) medium.
5. Thermal Stresses in Firn and Snow
Calculation of stresses
We shall now pass on to an analysis of thermal propagation and stresses in firn and snow. As pointed out in Section 3, the effective diffusivity of firn of density about 0.50 Mg $$ m-J is not significantly less than that of the high-density ice already treated, so the analysis of propagation of temperature waves can be carried over without modification. We can also retain the value quoted for the expansion coefficient a, since a connected structure of matter with air spaces expands at the same rate as a solid body of the same matter (it is just a matter of geometrical enlargement). Thermal strain-rates therefore remain the same. Since there is something of a contradiction in assuming that a firn layer could be subjected to the temperature wave of our Fossil Bluff model and yet remain firn, when in fact Fossil Bluff is a region of extensive summer melt, we shall also look at the effect of surface temperature waves corresponding to more realistic lower mean annual temperatures. We shall do this simply by reducing the term To in Equation (5) from -11°C to -21°C and to -31°C, but retaining the other parameters relating to the fluctuation of temperature about this mean
We shall not be able to follow the stress analysis of Section 4, since (a) the rheology of firm is not the same as that for ice, and (b) we shall find that we can no longer assume, as we could for high-density ice, that no tensile fracture occurs. To show this let us look at the flow and fracture behaviour of firn, and in particular its behaviour under tension.
Detailed study has been made of the densification of snow under compression, which can occur by a number of mechanisms, for instance: melting and refreezing; grain-packing by rearrangement; evaporation and condensation; grain-boundary diffusion of molecules; volume diffusion of molecules; plastic flow; and viscous flow (Anderson and Benson, 1963).
In the case of snow, densification proceeds primarily by the process of readjustment and closer packing of grains, while simultaneously there is the continuous process of "spheroidization" of grains, which occurs by processes such as evaporation and condensation, diffusion, and plastic flow. Grain packing continues until a critical close-packing density is approached, at a porosity of around 36-40%. After this point densification must proceed by mutual intrusion of particles (Benson, 1962).
Feldt and Ballard (1966) have summarized the laboratory and field results for the consolidation of snow and have provided a theoretical formula in good agreement with these results which treats snow as subject to Newtonian (viscous) flow characteristics and as having a macroscopic coefficient of viscosity calculable according to:
where p is the porosity of the snow, and b and $$ φ are empirically fitted constants. This form of relation is derived on the assumption that flow is dominated by the process of viscous flow of inter-particle bonds, and proceeds according to this until the critical density is reached. The critical porosity of 36-40% corresponds to density between 0.59 and 0.55 Mg m $$ K
The parameter N/V is strongly temperature dependent, so that snow becomes rapidly more viscous at low temperatures; experimental data from various sources are given by Feldt and Ballard (1966, p. 156), and are plotted here in logarithmic form against temperature (Fig 14). A straight line has been fitted to this logarithmic plot. The data show a large degree of scatter, owing to the many factors which can affect the viscosity of a snow sample of given density and temperature; its method of formation or preparation, its degree of age-hardening, its grain-size, and also the rate of loading used when testing it. The best value for b was found to be 1.78.
Fig. 14. Logarithmic plot of n/v against temperature, to determine the relationship between snow viscosity and temperature. Data from Feldt and Ballard {1966).
Literature on the tensile behaviour of snow is sparse. It has however been shown (Shinojima, 1967) that behaviour under tension is only slightly different from that under compression. He finds, admittedly in a different range of densities and stresses, that the coefficient of viscosity under tension is typically greater than that under compression by a factor of 1.61 but shows the same dependence on density and temperature. The greater viscosity under tension has been explained (Salm, 1975) as due to the fact that the mechanism of flow assisted by pressure melting under compression is no longer available. If we combine all these findings we can construct a general empirical formula for the effective macroscopic viscosity ηβ of firn under tension:
or
where Θ is the temperature in degrees Celsius. On this basis we can proceed directly to the stresses resulting from the strain-rates caused by temperature variations. In considering the possibility of fracture under tension, we shall be concerned not with the average stresses occurring in the firn cover, but with the maximum stress occurring during the year; for a crack to form we only require one occasion on which the critical fracture stress is exceeded. We calculate the stresses using the strain-rates used in Section 4 and the viscosity of Equation (20). Table I shows the maximum stress calculated in this way, and also a "typical" value of high stress occurring during the winter period, when most of the high stresses occur. By "typical" we mean here a value of stress which is exceeded on 50 d out of the 365 d during the year, and is therefore quite a common occurrence. These stress values are calculated for three different temperature regimes : those with mean annual temperatures of — 11°C-21°C and -31°C
The effect of temperature reduction is dramatic, to the extent that, if we assume no cracking occurs, then stresses at the surface attain absurdly high levels. If indeed no fracture did occur then we should expect these large stresses, alternately tensile and compressive, to play an influential role in the diagencsis of snow, for they are well in excess of the stresses due to overburden of snow; we could picture the surface cover of a glacier as undergoing a vigorous alternating stress, which must result in accelerating of grain rearrangement and more rapid firnification. Let us however consider the possibility that fracture occurs and relieves the stress.
Table I. Thermal stresses set up in firn, for three values of mean annual temperature. The division line shows the depth to which crack intiation is possible
The tensile strength of firn has been studied (Butkovich, 1956) and is found to be a function of density and temperature : the fracture stress is found to increase quite rapidly with increasing density, and to depend on temperature only by a factor of about 0.9% deg"1. For firn of density 0.50 Mg m-3 and temperature between — io°G and — 40°G we find the tensile strength to lie in the range 4.1 to 5.2 bar.
Table I shows stresses present as a function of depth, and we sec that for T0 — —11 °C tensile cracking can be initiated down to just over 0.5 m depth; for To = — 21°C it can be initiated down to 1 m depth; and for T0 = —31°C it can be initiated down to 3 m depth. Once initiated, we should expect the process of stress concentration at the tip of a crack to allow it to propagate well beyond the region in which the fracture stress is achieved (Parmerter, 1975) and its eventual depth will be governed by the same rules as apply to crevasses (Weertman, 1977). Indeed the thermal crack will become a crevasse, and this observation leads to a possible explanation of the initial stages of formation of a crevasse, for crevasses are often found to form in places where the tensile stress driving gross deformation does not exceed 4 bar, which is the minimum stress required to initiate a crevasse. The very high thermal stresses are enough to "nucleate" crevasses, which can then grow if the gross stress conditions are suitable.
6. Evidence
There is certainly direct evidence that the propagation of temperature waves through snow and ice can give rise to more than just gentle expansion and contraction: Swithinbank (1957, p. 25) remarks on narrow cracks in surface snow observed only during winter on the Maudheim Ice Shelf and suggests that they are the result of thermal contraction; seismological observation in Antarctica (Kaminuma and Takahasi, 1975) have shown that the propagation of the nocturnal cold wave gives rise to pronounced seismic events as the ice expands and contracts, especially when the temperature drops below —35°C. This is evidence that the formation of thermal cracks is likely, but the hypothesis of thermal straining and cracking also provides a satisfactory explanation of some otherwise unexplained experimental phenomena. Robin (1958, p. 48-50) and Kapitsa and Sorokhtin (1963) experienced unusually high seismic background noise when carrying out seismic shooting in regions of surface temperature below about —3O°C. If we consider our calculations of stresses to be properly significant only below about 1 m depth, owing to the considerations mentioned in Section 3 and the fact that the surface snow is largely recent and unconsolidated, then our Table I provides an explanation of why significant fracture occurs only at temperatures below about — 3o°C. It has also been found (Bentley, 1964, p. 344) that the amount of seismic noise at "Little America" station and at "Byrd" station shows seasonal variations, and is at a minimum during the early summer. The mean annual temperatures at these stations are — 22°C and — 28°C respectively. It is possible that thermal cracks close up and weld themselves together during the summer, and only form again as winter approaches and stresses again reach the fracture strength. On the high plateau, at mean annual temperatures of about -6o°C, the fracture strength is reached throughout the year and seismic shooting noise is present always.
We shall now try to say something about the width of a purely thermal crack and the character of the strain-rate field surrounding it. Let us consider a field of thermal cracks at spacing λ (Fig. 15). They will open and close as temperature fluctuates, but the width of opening depends on the degree to which ice between cracks is free to expand and contract. At the very surface, immediately on the boundary of a crack, the ice is free to expand and contract in response to fluctuations in air temperature, but any disturbance is attenuated with depth owing to the conduction properties already analysed and is also attenuated with distance inwards from a crack according to Saint Venant's principle concerning edge effects. Depth attenuation is frequency dependent, and occurs according to exp [—j(u>/2A)']. We can therefore define a "skin depth" j0, in analogy with electromagnetic theory, as yn = (2k/w)1 representing the depth at which temperature effects are attenuated by a factor 1/c. We shall use this skin depth as a characteristic distance for the attenuation of thermal effects in the y direction, and shall make the further assumption that it also characterizes the attenuation with distance inwards from a crack, in the x direction. This is a way of expressing the principle that edge effects generally die away beyond the order of a "few thicknesses" from a free boundary. This means that we can approximately represent the strain-rate field as a function of x and y by
where {ixx)a is the strain-rate which a free ice mass would undergo if subjected throughout to the surface air temperature fluctuation, and
is simply the distance away from a crack.
Fig. 15. Diagram of the strain-rate field surrounding cracks undergoing thermal strain; appreciable strain occurs within a region of order j0 surrounding a crack. A wire strainmeter with its end units in the vkinity of a crack will register an anomalous apparent strain-rate.
We can now calculate the rates of opening of cracks: we take the origin of coordinates as a stationary point by symmetry considerations, and then integrate strain-rate with respect to x to find the rate of movement ve of a crack edge. This gives, for the surface,y = 0,
and provided that A > y0, i.e. that the crack spacing is very much greater than the skin depth, we can write
meaning that the velocity of opening and closure of a crack is:
For the final maximum width W of a thermal crack we can write
where (eXx)max is the maximum strain of a free ice mass. This reduces to:
To make quantitative estimates we shall need values for skin depth. Many different frequencies of thermal fluctuation are present, with many different associated skin depths, but for the purposes of round numbers we shall look first at the case of a surface temperature oscillation of amplitude 7 deg and of period about one month, for which the associated skin depth is about one metre. We find that:
and
and therefore
The maximum width of opening, in the absence of enhanced opening due to gross deformation, is:
These figures refer to a temperature wave with skin depth 1 m, so may be taken to give the order of magnitude of crack parameters in the region of the top metre and over a period of about one month ; shorter-period waves give rise to more rapid, but smaller, opening and closing, and longer period waves give rise to slower, but larger, opening and closing. We expect the true maximum width of a crack to be determined by the cycle of longest period (which has the greatest skin depth and hence the greatest horizontal extent). An annual cycle of amplitude 11 deg, and skin depth 3.3 m, gives rise to a maximum crack width of 7.7 mm. A crack of width the order of a few millimetres is easily wide enough to present an effective reflector to seismic waves, for the amplitude of firn oscillation due to the passage of a seismic disturbance is of a very much smaller magnitude, and is insufficient to bridge the air gap.
Thermal straining is also of sufficient magnitude to provide a possible explanation of some anomalous effects encountered in work with recording wire strainmeters (Evans and others, 1978). In this work, strainmeters with gauge length 5 m were set up on the Barnes Ice Cap with a view to evaluating the consistency of strainmeter results with results obtained using more conventional stake networks over longer gauge lengths. It was found that strainmeter readings were subject to an unexplained fluctuation with a period of some 11 d. Strainmeters were set up on solid ice after clearing 1 m of surface snow, and the temperature of the strainmeter was monitored. Insufficient temperature data were recorded at the time to be able to establish a definite connection between temperature variations and strain-rate fluctuation, but it is extremely likely that it was caused by some effect connected with weather conditions in the area, since the fluctuation shows up practically simultaneously on the traces of three independent and quite widely separated strainmeters (some 10 km apart). It cannot therefore be simply a local strain-rate anomaly, and the hypothesis that the fluctuations are due to the passage of a kinematic wave leads to unreasonably large wave velocities. It is remarked in the paper that "cracks of up to 1 cm width, and exhibiting progressive opening, were observed. The snow pack was split above most cracks, indicating opening during the winter. . . orientation tended to be variable”. The facts of irregular orientation and opening during the cold of the winter are quite strong indicators that the cracks are of thermal origin, although the fact that they occur in high-density ice contradicts our earlier assertion (Section 4) that thermal stresses cannot achieve the fracture strength of high-density ice. However, the origin of the cracks is irrelevant to their subsequent flexure under thermal stress. The spacing of the cracks was of order 5 m and the strainmeters were always set up so that although they were often close to cracks their line of action never actually crossed any. Gross strain-rates observed were of the order of 10-10s-' and fluctuating strain-rate was of amplitude about 2X101ls -1 about this value, with an apparent period of about 11 d, although all these statistics are based on a data set covering only 21 d.
Consider a strainmeter set up as in Figure 15, with the two end units set up at distance & away from a crack. By integration of kxx between the stationary origin and the strainmeter end units we find that the velocity vs of an end unit is given by
or
for
The relative velocity of the two end units is therefore
and the apparent strain-rate registered, calculated for the entire gauge length of λ--2δ, is
The distance δ is unfortunately not known, and is also a critical parameter in the expression for kxx. However, since the crack spacing on the Barnes Ice Cap was of the order of 5 m and this is also the gauge length of strainmeter being used, it is fair to say that strainmeter end units must often have been in the "vicinity" of cracks—that is within one skin depth of them. For a wave of period 11 d the skin depth is 0.58 m, and we shall carry out the calculation for a wave of amplitude 3 deg. corresponding to that recorded on the Barnes Ice Cap. So, we take δ as 0.58 m, and take the gauge length to be 5 m. Then
may be taken as a rough estimate of the apparent strain-rate due to thermal effects. It is seen that this is easily of sufficient magnitude to account for the fluctuation observed, and is indeed large enough to account for all the strain-rate observed, including the gross strain about which the fluctuation occurred. It is unlikely for the case of the Barnes Ice Cap that the gross strain-rate is connected with thermal activity, since conventional, long gauge-length strain rosettes were found to give fairly acceptable agreement with strainmeter results. It is clear nonetheless that thermal straining in the vicinity of cracks can be of at least the same order as gross dynamic strain, and steps should be taken to avoid it by use of long gauge-length strain networks or strainmeters set up well away from any crack, or to eliminate it by taking enough measurements (of temperature, cross-correlated strainmeters, crack location and activity, etc.) that temperature effects are known and therefore subtractible.
7. Conclusions
Let us summarize our findings:
High density ice (Tu = — 11°C
(i) Thermal stresses are appreciable within the top 3 m and are generally in the range0.5 to 1.5 bar.
(ii) Fracture cannot occur unless assisted by stress concentrations.
(iii) Conditions are favourable for the formation of surface rumples of wavelength about 10 m and initial amplitude less than 20 cm. This is therefore not an explanation for the melt pools described.
Firn and snow regions (T0 below about —20°C)
(i) Thermal stresses are dependent on temperature and often exceed 4 bar in the top 3 m.
(ii) Fracture can certainly occur, and propagate itself.
(iii) Thermal stresses and strains provide explanations for:
(a) observations of surface cracking during cold weather ;
(b) seismic shooting noise in cold polar regions;
(c) initiation of crevasses in regions of relatively low overall stress;
(d) anomalous strain-rate measurements near cracks.
8. Acknowledgements
I should like to thank Dr C. S. M. Doake for useful discussions, and Messrs J. F. Bishop, J. L. W. Walton, and R. I. Kressman for carrying out the computation involved. I am grateful to Dr G. de Q. Robin for drawing my attention to problems encountered in seismic shooting, and to Dr D. J. Goodman for telling me of his results on the Barnes Ice Gap.
