Thermal behaviour of glacier and laboratory ice

Abstract

Water is present in glaciers in the form of veins at the three-grain junctions. This water remains unfrozen even many degrees below the normal freezing point, mainly because it contains much of the soluble impurity content of a glacier, but also because of the microscopic curvature of the ice–water interfaces. As the temperature is lowered and the veins shrink, the concentration of impurities in them increases, and the curvature effect also increases. The predicted relation between vein size and temperature has now been verified by laboratory experiments.

Because of the latent heat of the vein water, the ice behaves macroscopically as a continuum with an anomalous specific heat capacity that depends strongly on temperature. From this point of view, a polythermal glacier is a single medium with continuously varying properties, rather than consisting of distinct cold and temperate phases with sharp boundaries between them. The paper sets up differential equations for heat diffusion in such a continuum. To explain the local uniformity of the vein system seen under the microscope, it is found necessary to include the effect of diffusion of solutes along the veins.

Solutions are presented for a model in which two semi-infinite slabs, initially having different temperatures, impurity concentrations and vein sizes, are instantaneously brought into contact. In this way, transition thicknesses between cold and temperate ice are estimated, and also the velocities of various kinds of waves that are generated from the original discontinuity at the interface.

1. Introduction

One approach to the problem of a polythermal glacier (Fowler and Larson, 1978; Hutter, 1982; Fowler, 1984; Hutter and others, 1988) is to model it as consisting of two distinct continuous phases, namely, temperate and cold (or polar) ice. Temperate ice is at the melting point and contains moisture, while cold ice is below the melting point and is dry. Appropriate differential equations are formulated for the two different phases, and jump conditions are defined for the interface between them.

We adopt here the different point of view that polycrystalline ice can be modelled as a single medium, but with smoothly varying thermal properties. Of course, this is not necessarily in conflict with the two-phase model, which may still be appropriate if the transition layer between the cold and temperate phases turns out to be sufficiently thin. The view of ice as a single medium with smoothly varying properties is based on observations in the laboratory of the behaviour of single veins of water in polycrystallinc ice at the three-grain junctions (Mader, 1990, in press a, b). At temperatures within the range 0° to about – 1°C, change of temperature is observed to be accompanied by a change in the size of the veins. This results from two effects (Lliboutry, 1971; Harrison, 1972; Nye and Frank, 1973; Raymond and Harrison, 1975; Harrison and Raymond, 1976).

The first arises from the impurities that are always present dissolved in the vein water. Because the ice phase tends to reject impurities, the shrinking of a vein by freezing concentrates the salt solution it contains, and so lowers the equilibrium temperature. The second effect arises because the vein walls are highly curved. The narrower the vein, the greater is the curvature and, because the vein walls are concave inwards Fig. 1, the narrower the vein the lower is its equilibrium temperature. Thus, both effects work in the same direction.

Fig. 1. Cross-section of a water vein at a three-grain junction in polycrystallinc ice.

Lowering the temperature does not freeze all the vein water; rather it progressively reduces its amount. The latent heat of the vein water means that, viewed as a continuum, the ice appears to have an anomalous specific heat capacity, which becomes large as the veins become large. The anomalous specific heat capacity gives rise to an anomalous thermal diffusivity, which determines how long an inhomogeneity of temperature of given spatial scale can persist.

The concentration of impurities within the veins changes not only because of changes in the vein size, but also by virtue of diffusion of impurities along the vein network, in response to gradients of concentration. Although this diffusion process is slower by a factor of 10³ than the unmodified thermal diffusion, it has to be considered if we are to account for the local uniformity in vein size that is observed, as we shall see in detail later.

How all these processes interact together is not at first obvious. Therefore, the primary purpose of the paper is to study the continuum thermal behaviour implied by the observed vein structure within the ice. To this end, we first formulate the appropriate partial differential equations for the ice viewed as a continuum, and then use them to study how uniformity (of temperature, concentration and vein size) is reached from given starting conditions. In particular, we calculate the time-scales and length-scales that are involved. There are some similarities with the solidification problem in alloys (Worster, 1986).

[Since the submission of this paper the same continuum model has been used to study the slow coarsening of the veins that takes place in ice near its melting point when stored in the laboratory (Nye, 1991b).]

Before applying the results directly to glaciers, several further points would need to be considered. For example, glacier ice may contain some proportion of water that is not at the veins, in the form of lenses at grain boundaries and in irregular shapes (Raymond and Harrison, 1975). This will contribute to the specific heat anomaly; however, the fact that the veins are the sites where the water is in thermodynamic equilibrium suggests that they may often contain most of it.

Another effect to consider (in laboratory experiments as well as in glaciers) is the movement of water along the veins. This could occur in at least three ways. (1) In a glacier, water, being the heavier phase, will tend to flow downwards under gravity. (2) The change of specific volume of ice on melting or freezing of veins can cause pressure gradients and consequent flow of water. (3) Vein water will also flow under the influence of other pressure gradients, generated, for example, by the hydraulic system in a glacier or by ice movement.

The relative importance of these processes will depend on the precise physical conditions. Process (l), with the water conveying heat, is usually considered to be the principal mechanism by which a glacier becomes “temperate” despite being cold in its higher reaches. This effect will depend on the direction of the temperature gradient with respect to the vertical. As another example, process (1) would not operate in an experiment where a block of polycrystalline ice in the laboratory is held totally immersed under water, so that the vein water is simply in hydrostatic equilibrium. However, such a block will be subject to process (2) if its temperature is lowered so that the veins begin to freeze; the expansion would be expected to drive water along them out to the surfaces, and this process has been observed (Mader, 1990, in press b). Process (2) may not occur so readily within a much larger body of ice such as a glacier, where there will usually be no nearby free surfaces to which the expansion water may flow. If water flows, there will also be dispersion of the solute in addition to diffusion (Bear, 1979). Thus, so far as water flow along the veins is concerned, it would appear that each physical application has to be considered on its merits. It is for this reason that we confine ourselves here to the well-defined problem where water flow along the veins is not allowed.

Our model also neglects the effect of bulk strain rate, although Robin (1955) has shown that this can be an essential element in maintaining a steady distribution of temperature within an ice sheet. For all these reasons, one should be cautious in drawing firm conclusions about glaciers from this simple model. Its purpose is to deal precisely with one aspect of a complex problem.

2. Continuum Model

2.1. The Relation between Temperature and Vein Size

Consider a single vein containing a concentration c (mass per unit volume of liquid) of soluble impurities. Let the radius of curvature of its cylindrical faces be r _v (Fig. 1). In equilibrium, the vein and the surrounding ice will be at the melting point, which is depressed by

(1)

where v_i, is the specific volume of ice, s_w and s _i are the specific entropies of water and ice, γ is the ice–water surface energy and Κ is a positive constant that depends on the impurity composition (Schwerdtfeger, 1963). It is convenient to define an equivalent vein radius r such that a vein of circular section and radius r would have the same cross-sectional area: πr ² = αr _v ².

The value of α depends on the dihedral angles at the vein edges, that is, the angles ϕ of the curvilinear triangular cross-section of the vein. Measurements by Walford and others (1987) of the focal lengths of water lenses at grain boundaries gave ϕ = 33.6 ± 0.7°. This deduction used simple lens theory based on geometrical optics. However, more accurate diffraction theory (Nye, 1991a) has revealed a systematic error in the method arising from diffraction, and new measurements of the diffraction pattern (Walford and Nye, 1991) have given ϕ = 25.0 ± 1.0° for a supposedly typical grain boundary. At the same time, there is evidence (Morris, 1972; Mader, 1990, in press a) that ϕ can vary considerably with the relative crystal orientation at the grain boundary. We shall provisionally take for calculations ϕ = 25.0°, which corresponds by the formula (1) in Nye and Mae (1972) to α = 0.1007.

Equation (1) may now be written, using r, as

(2)

where

With ν_i = 1.09m³Mg⁻¹, γ = 34mJm⁻², s_w – s _i = 1.22 kJ kg⁻¹ deg⁻¹, we find A = 5.44 × 10⁻⁹m deg. Equation (2) expresses the depression of the melting point u in terms of vein size r and impurity concentration c

2.2. Conservation of Energy

We consider a polycrystalline sample of ice containing a vein density λ, in units of vein length per unit volume. In order of magnitude, λ will be related to the grain-size b by λ 3b⁻². The veins have equivalent radius r and the temperature depression is u. If r and u now change by dr and du, we write the heat required per unit volume of the composite as

(3)

the first term arising from the latent heat of the vein water (L per unit mass) and the second from the specific heat capacity of the ice (σ per unit mass). Consistent with our assumption of no water flow along the veins, we neglect the volume change on melting and use a common density ρ throughout. In expression (3) we are taking the relative volume of the veins to be small, and we are neglecting the contribution from the specific heat capacity of the vein water, in comparison with its latent heat. The former assumption is well satisfied in practical cases, and the latter is valid provided the total temperature change considered is much less than L/σ = 160 deg, which is adequate for our purposes.

Expression (3) leads to the energy equation

where k is the thermal conductivity of ice, essentially unchanged by the presence of the veins, or

(4)

where D_th = k/ρσ, the thermal diffusivity of ice. This is the standard heat-diffusion equation, but with the addition of the second term, which is due to the veins. To see its physical significance, consider the case when no heat flows in. The righthand side is then zero and the two terms on the left balance. Thus, to increase the radius by melting requires heat, which leads to a fall of temperature (an increase in u).

2.3. Conservation of Mass of Solute

Let us first imagine for simplicity that all the veins are parallel and run in the x-direction. The flux of solute along a vein by virtue of the concentration gradient is –D _c∂c/∂x, where D_c is the diffusivity of salt in water. This is per unit area of vein. The volume of vein per unit volume of continuum is πr ²λ. So the flux per unit area of continuum is –D _cπr ² λ∂c/∂x. Since the mass of solute per unit volume of continuum is πr ²λc, this leads to the conservation equation

(5)

We generalize this for veins distributed is otropically to

(6)

where the factor 1/3 arises because individual veins are more or less favourably oriented with respect to grad c. To justify the factor, first note that the flow of solute, mass per unit time, along a vein inclined at an angle θ to grad c is proportional to cos θ. If the density of veins at angles between θ and θ + dθ (length per unit volume of continuum) is dλ, the number of such veins crossing unit area normal to grade is cosθdλ. Hence, the flow across unit area of continuum due to this class of veins is proportional to cos²θd λ. Since the factor cos²θ averages to 1/3, the result follows.

The implicit assumption in Equation (5) and (6) is that, when a vein freezes, all the impurities remain in the liquid phase; no salt from the vein is trapped, for example, in the three grain boundaries as they extend. This is an essential feature of the present formulation. It was originally adopted as a plausible working hypothesis, but the experimental work of Mader (in press b), which was specifically undertaken to test it, has now given it a firm empirical basis.

The three Equations (2), (4) and (6) determine the space and time behaviour of the three fields u, r and c.

2.4. The Anomalous Specific Heat Capacity

The idea of an anomalous specific heat capacity is appropriate when solute is not allowed to diffuse along the veins; for, if it were allowed, the input of heat would not be the only factor that determines temperature change. Therefore, to see how it is included in the equations we put D_c = 0.

Then Equation (6) gives

, and, if we denote by m the mass of solute per unit length of vein (m = constant), we have

(7)

Equation (2) is now

(8)

where the constant Β = Κm/π.

Since u is now a function only of r, and vice versa, the second term in Equation (4) involving ∂(r ²)/∂t can be expressed in terms of ∂u/∂t like the first one. Thus,

and Equation (4) can then be written simply as

with an effective diffusivity D _eff which is in turn related inversely to an effective specific heat capacity σ_eff. Thus,

(9)

or

(10)

The last term is the anomaly.

Equations (8) and (10), with r as a parameter, determine the effective specific heat capacity as a function of temperature depression u. As u → 0, r, → ∞ and the anomaly becomes infinite. For a given impurity content m and small enough vein radius r, the first term in Equation (8), arising from curvature, becomes negligible compared with the second, arising from impurities. In that case, putting A = 0, we can eliminate r between Equations (8) and (10), and obtain explicitly

(11)

where

(12)

in agreement with Harrison (1972, p. 26). u_t is a transition temperature depression. When u falls to u_t, the specific heat capacity becomes twice its normal value. This simplified formula is valid when u ≫ A ²/B or u ≫ πA ²/Km.

The behaviour of the effective diffusivity D_eff as a function of temperature depression u is sketched in Figure 2. It falls to zero as u → 0. Therefore, there comes a point where it falls below (1/3)D_c, which, from Equation (6), is the diffusivity that governs the movement of impurities along the veins. For lower values of u, the notion of the anomalous specific heat capacity as governing the diffusion of heat breaks down; it is in fact D_c rather than D_eff that then controls the approach to equilibrium, as we shall show by numerical solutions in section 6. We shall call such ice “deeply temperate”.

Fig. 2. Showing diagrammaticaLly the dependence of the effective thermal diffusivity D eff on t emperature. The temperature depression ut below the normal melting point (the origin) marks the cold–temperate (C–T) boundary. Where D_eff < D_c/3 the ice is deeply temperate (D.T.). Note the highly non-linear scale of diffusivity, which is in units of Dth.

2.5. Numerical Examples

Throughout this paper we shall use two numerical examples, chosen to typify laboratory ice and glacier ice. In example (1), a vein in a laboratory specimen of ice, grown from singly distilled water and purified by an additional factor of about ten during freezing (Mader, in press a,b), was observed to have r = 6 µm at a temperature depression u = 0.10 deg. In example (2), from Raymond and Harrison (1975), a specimen of fine-grained glacier ice was inferred to have, when in the glacier, r = 14 µm at at a temperature depression of 0.006 deg (point P ₃ in their Figure 9). (This docs not include the temperature depression due to pressure.)

So the curvature term in Equation (8) is A/r = 9.1 × 10^–4 deg for example (1) and 3.9 × 10⁻⁴ deg for example (2), which we neglect in comparison with the measured u. Then B = ur ² = 3.6 × 10^–12 deg m² for (1), and 1.2 × 10⁻¹²deg m² for (2). Taking λ = 5 × 10⁴ m⁻² for (1) and 1 × 10⁶ m⁻² for (2) and L/σ = 160 deg, gives u_t = 0.0093 deg for (1) and 0.024 deg for (2).

If these specimens were heated up so that the veins grew to a size where the curvature term in Equation (8) equalled the impurity term, the temperature depressions would be u = 2A ²/B = 1.6 × 10^–5 deg for example (1) and 5.0 × 10⁻⁵ deg for example (2). The vein radius would be 660µm for (1) and 220µm for (2). Thus, for these specimens, the simplified formula (11), which neglects curvature, is applicable at least up to temperatures where u ≈ u _t. u _t, which may be regarded as marking the transition temperature from cold to temperate, corresponds to a vein radius r = 20 µm for (1) and 7.0 µm for (2). Veins even as small as this soak up enough heat to generate a 100% anomaly in the specific heat capacity.

The temperature at which D_eff falls to (1/3)D_C (Fig. 2) can be estimated from Equations (8) and (9). For the fine-grained glacier ice of example (2), we find u = 5.4 × 10^–4 deg, which corresponds to a vein size r = 52 µm. For higher temperatures and larger veins than this the ice is deeply temperate.

3. Differential Equations

We shall be concerned with one-dimensional solutions and therefore we write the governing Equations (2), (4) and (6) as

(2)

(4’)

(6’)

the dependent variables being u, r and c.

To express them in dimensionless form, take the unit of temperature depression as L/σ = 160 deg, the unit of length as , which is of the order of the grain-size, and the unit of time as (λD _th)^–1. It is now more convenient to work in terms of the area rather than the radius of the vein cross-section. Denoting dimensionless temperature depression, distance, time and area by U, Χ, Τ and S, we have

S is then the fractional vein volume. Further, we define dimensionless concentration C by

Then the equations for the fields U(X,T), S(X,T) and C(X, T) reduce to

(13)

(14)

(15)

where just two dimensionless constants remain, namely,

.

With D _c = 1.3 × 10^–9 m² s^–1 and D _th = 1.2 × 10^–6 m² s^–1, their values are

(A further scaling which has the effect of removing a altogether could be accomplished by defining new variables obtained by multiplying U, S and C all by a^−2/3. However, we prefer to work with dimensionless variables whose physical meanings are more apparent.)

4. The Two-Slab Model

A simple model, designed to show the thermal behaviour, and particularly the time and space scales, implied by Equations (13), (14) and (15), is the following. Two semi-infinite slabs of ice, with a common interface X = 0, each of them initially uniform but having different temperatures, vein sizes and impurity concentrations, are suddenly brought together at Τ = 0 and allowed to equilibrate. This problem allows a similarity solution.

Define the new variable

(16)

and hypothesize that, for our initial value problem, U, S and C are functions not of X and Τ separately, but solely of Z. We shall write them as U(Z), S(Z) and C(Z).Then the partial differential equations become the ordinary ones

(17)

(18)

(19)

where primes denote differentiations with respect to Z.

The boundary conditions in (X,T) can also be expressed in terms of Z, namely, that U(∞), U(-∞), C(∞) and C(—∞) should have prescribed values. These boundary conditions are appropriate to the two-slab problem because, at Τ = 0 and for all X > 0, Ζ = ∞, and at Τ = 0 and for all X < 0, Ζ =-∞. Thus, U(∞), U(-∞), C(∞) and C(-∞) are the initial (T = 0) uniform temperature depressions and concentrations in the two semi-infinite slabs. Figure 3, which is a key to interpreting all the graphs that follow, illustrates this. Having fixed these values, the initial values of S in the two slabs are determined by Equation (17).

Fig. 3. A sketch graph to illustrate the trans formation Ζ = X/T^½. The curve shown for U(Z) may be read in two ways. (1) For fixed Τ it shows the spatial distribution along X of temperature depression U. As Τ increases, the curue becomes uniformly stretched out, by a factor T^½. (2) For fixed position X the curve shouis how U changes with time. For X > 0, time runs from Ζ = ∞ to Ζ = 0, and for Ζ < 0 time runs from Ζ =-∞ to Ζ = 0. Thus, for X > 0, U relaxes from. U(∞) to U(0), and, for X < 0, U relaxes from U(—∞) to U(0).

It is useful to take, as a subsidiary dependent variable, the mass m of solute per unit length of vein, which was defined in Equation (7) but there held constant. If we define the dimensionless mass per unit length by

(20)

Equation (7) is simply

(21)

In this model, M can only be changed by diffusion of solute. Any pair of the variables U, S, C, M (or u, r, c, m) suffices to determine the state of the ice.

It is an immediate and interesting consequence of the transformation from (Χ, Τ) to Ζ that, after the first instant (T = 0), the common interface X = 0 remains fixed in temperature throughout the diffusion process.

The limit Τ → ∞ merits further comment. Τ = ∞ corresponds to Ζ = 0 for all finite X, but not necessarily for X = ±∞. There are two limits involved here, Τ → ∞ and also the fact that the slabs are semi-infinite, and the order in which they are taken is important. If the thickness of the slabs is allowed to be infinite first, the temperature distribution always has the form calculated in our model; for fixed T, however large, there are sufficiently distant parts of the slabs that retain their original temperatures. If, on the other hand, we started with slabs of finite thickness and let Τ → ∞, the development of the temperature distribution would be different, because ultimately the distant faces would have an influence. If the distant faces were insulated, the uniform temperature ultimately attained would depend on the relative thicknesses of the slabs, and would not in general be U(0). This would continue to be so even if the thicknesses were then allowed to become infinite. However, for finite slabs, the temperature depression at the interface for Τ small (but not zero) would be U(0), because in this case they would act as if they were semi-infinite. In practice, it need hardly be said, the sizes of blocks of ice are finite, and so are their lifetimes. For a glacier, the space and time scales make the model of semi-infinite slabs an appropriate one. In the laboratory, on the other hand, experiments could be devised to suit either model.

5. Behaviour of the Equations in the Limits D_c → 0, A → 0

At temperature depressions greater than 0.1 deg, the contribution to temperature depression made by curvature is less than a few percent, and so one might be inclined to neglect it for many purposes, retaining only the contribution from salinity. One would put A = 0 in Equation (2), which is equivalent to setting a = 0 in the dimensionless Equation (17). Also, because the diffusivity for salt in water D_c is some 10³ times smaller than the diffusivity for heat in ice, one might be inclined to put D_c = 0, which is equivalent to setting d = 0 in Equation (19).

If one tries to simplify and solve the two-slab problem in this way, there is no difficulty, except that S and M arc, in general, discontinuous at Ζ = 0. The temperature depression U is always continuous, and so is the concentration C, because in this case U = C. This means physically that, even at t = ∞, the vein size is discontinuous across t = 0, the concentration is continuous, but the mass per unit length is discontinuous. Thus, putting a = 0, d = 0 leads to no mathematical inconsistency.

Suppose now that we try to relax this approximation by putting a ≠ 0, d = 0 that is, allow for a curvature effect but continue to neglect diffusion of solute. There will now be, in general, discontinuities at Ζ = 0 not only in S and M but also in C. This again is quite consistent (and we show numerical results in section 6).

It might be supposed that all these discontinuities would be removed if d ≠ 0, to allow diffusion of solute, and indeed they are. However, it is perhaps surprising that they are not removed if we deal with the case a = 0, d ≠ 0. A physical argument is best for explaining why this is so. In Figure 4, where a = 0, d ≠ 0, a vein has discontinuities in cross-sectunion and in mass of solute per unit length, but the temperature depression and the concentration are both uniform. All the equations are satisfied and there is equilibrium. The non-zero value of d does not suffice to smooth the discontinuity in the vein because the concentration, which alone controls the temperature (a being zero), is uniform. However, if we now suddenly allow a ≠ 0, we generate a temperature difference, heat will flow and the vein becomes uniform. Thus, in summary, the case a ≠ 0, d = 0 and the case a = 0, d ≠ 0 both allow discontinuities. Therefore, both diffusion of impurities and the effect of vein-wall curvature, small though it may be, are essential for attaining the local uniformity of the veins that is seen under the microscope.

Fig. 4. Αn idealized longitudinal section of a vein in equilibrium in a model where the curvature of the vein walls has no effect on melting point (a = 0), and where the diffusivity of solute is non-zero (d ≠ 0). There can be discontinuities in cross-section S and in mass M of solute per unit vein length, but the temperature depression U and concentration C are both uniform.

6. Solutions of the Equations

We note from the definition Ζ = X/(T)^½ that, for all finite Χ, Ζ = 0 corresponds to Τ = ∞. Thus, the final uniform values of U and C are given by [U(0) and C(0). To integrate Equations (17), (18) and (19) numerically, the procedure adopted was to set values for U(0), C(0) and also for U’(0) and C’(0), and to integrate outwards for positive and negative Ζ until asymptotic values U(∞), U(—U), C(U) and C(-U) were reached.

In practice, it was hard to find a method of numerical integration that was always stable. Because of the small value of d (3.6 × 10⁻⁴), there can be fine detail in the functions for comparably small values of Z. Physically, this reflects the two widely different diffusivities involved. In many cases the thermal diffusivity D_th dominates when |Z| ≈ 1, while D_c is effective when |Z| ≈ 10^–3. The numerical schemes adopted are discussed in Appendix A.

To illustrate the role of D_c, computations were also made (Appendix B) with the same asymptotic values as before, but setting D_c = 0 or equivalently d = 0.

If we choose C′(0) = U′(0), it follows from the governing Equations (17), (18) and (19) that S′, S″, C″, U″, M″ are all zero at Ζ = 0. Thus, this choice will correspond to solutions that are rather smooth at Ζ = 0. By contrast, if we chooseC′ (0) ≠ U′ (0), we can expect more detail at Ζ = 0, which will correspond to slow diffusion processes. This distinction makes it convenient to divide the example solutions into two groups: (A) with C′(0) = U′(0) and (B) with C′ (0) ≠ U′ (0). Numerical values for the five examples are collected in Table 1.

Table 1. Temperature depressions, vein sizes and impurity contents

(A) Solutions withC′(0) = U′(0)

We choose the value of M(0) (mass of solute per unit vein length at Τ = ∞) and the constant a to be those corresponding to the impurity content and grain-size of example (2) of section 2.5, which was fine-grained glacier ice. Thus M(0) = 2.30 × 10⁻⁸, a = 6.0 × 10⁻⁸. Then, from Equations (21) and (17)

Choosing U(0) now fixes S(0) and hence C(0). Thus, since C′(0) = U′(0), we are free now to choose as boundary conditions only U(0) and U′(0).

Let us agree to call ice temperate or cold, as suggested by Harrison (1972), according to whether its temperature depression is below or above the transition temperature depression u_t defined by Equation (12). In other words, ice is temperate if it is above the temperature at which the direct and phase-change contributions to the effective heat capacity are equal. By Equation (12), u_t depends on the mass m of solute per unit length of vein. If we denote by U _t the dimensionless temperature depression corresponding to u _t, so that

Equation (12) is simply, in view of Equation (20),

(22)

Thus, in terms of dimensionless variables, the definition is:

In broad terms there are three interesting cases according to the initial temperatures of the two slabs, and these translate into conditions on U(0) and U′(0).

• (i) Cold: cold. Here U(0) is large, so that the contribution from curvature is insignificant, and, because of the smoothness at Ζ = 0 (Τ = ∞), it turns out that D_c has negligible effect. We are in the region where it is a valid approximation to put a = 0, d = 0, so that U(Z) = C(Z). Equations (17), (18) and (19) then reduce to the single equation

  

  This is what would be obtained from the ordinary heat-diffusion equation, but with an anomalous diffusivity as discussed in section 2.4.

  If, in addition, U>>U_t, the anomaly term is negligible and we recover the standard diffusion equation with constant diffusivity,

  

  which has solution

  

  This is symmetric in the sense that the final temperature depression is the mean of the two initial temperature depressions.

  A typical solution of the full Equations (17), (18) and (19), without approximation, for this cold-cold transition is shown in Fig.5. The transition width, suitably defined as

  

  has the value 3.3. In terms of physical variables this is a width , which is 200 m at t = 100 year and 640 m at t = 1000 year. For a laboratory-scale experiment we find ∆x = 22 cm at t = 1 h.

• (ii) Deeply temperate: deeply temperate. We have chosen for this case the opposite extreme where the initial temperatures of both slabs are in the deeply temperate range; the effective thermal D _eff is smaller than Dc/3 so that the latter dominates. In Figure 6, which illustrates this, separate full curves are shown for U, C and a/S ^½. According to Equation (17), the dimensionless quantities C and a/S ^½ may be regarded as the contributions to U made by concentration and curvature, respectively.

  The transition width for U is now much smaller than in the cold-cold example (note the different Ζ scale), being ΔZ = 0.068, which is 4.2 m at t = 100 year, 13 m at t = 1000 year and 4.5 mm at t = 1 h. The smallness of ΔΖ is due to the small value of D _c; the ice takes longer to attain equilibrium and the U curve leaves the asymptotes only at small values of |Z|. A notable and unexpected feature is that the vein cross-section S is almost constant, with the temperature closely following the concentration. (Putting S′ ≡ 0 in the governing Equations (17), (18) and (19) leads to inconsistency. In fact, a small nonzero value of S′, barely visible on the graph of a/S^½, is essential for satisfying the equations.)

  If we had assumed D _c = 0, the curves for U, C and a/S^½ would have been those shown by the broken lines. The transition width is even smaller, a consequence of the very small value of D _eff. The broken curves are discontinuous at Ζ = 0. With D_c = 0 there would be discontinuities in C and S at the interface, even at Τ = ∞, as noted in section 5. The temperature would ultimately be uniform and continuous, but there would be no physical process to make the concentration and vein size continuous. In this model (which ignores any water flow through the veins), diffusion of impurities is essential if the vein size is to vary smoothly without discontinuity. In fact, in this example, as noted above, the whole process is dominated by D _c because both slabs remain deeply temperate.

• (iii) Deeply temperate: cold. This is the most interesting of these three cases and is illustrated in Figure 7a,b and c. We have chosen to have the righthand slab initially at −1.04°C (U(∞) = 6.5 × 10⁻³). The lefthand slab is deeply temperate with initial temperature depression U(—∞) = 1.03 × 10⁻⁶. This temperature depression corresponds, with our standard value of impurity level Μ, to a vein radius of r = 97 µm. The target values of (U(∞), U′(-∞) having been set, the values of U(0) and U′ (0) necessary for starting the integration were adjusted accordingly (U(0) = 7.00 × 10⁻⁴, U′ (0) = 3.40 × 10⁻³).

  The great range of scales makes it necessary to display the results on three graphs at successively higher magnifications. Figure 7a shows that on the cold side the behaviour is like one-half of Figure 5, with a Ζ scale of order 1. With our standard value of M(0), the transition temperature depression U _t, calculated from Equation (22) as , is 4.6 times less than U(0). So the common boundary in this example is cold. Figure 7b shows that the temperate-cold boundary U = U_t lies at Ζ = Zt = —0.16; that is, its equation of motion is , or, in terms of physical variables,

  

  The boundary moves from x = 0 to x = −9.8 m after 100 year and to x = —31m after 1000 year. In 1 h the movement is 11 mm.

  The curve of U(Z) in Figure 7b is remarkably linear, although, if one thinks in terms of a Taylor expansion about Ζ = 0, this is perhaps not so surprising in view of the small values of |Z|. However, U cannot become negative. This conflict results in a very sharp change in slope as U approaches zero, of which the detail is shown in Figure 7c.

  Physically, the curves of Figure 7a, b and c may be interpreted as showing that the cold side succeeds eventually in making the whole composite slab cold. It does this by moving the temperate-cold boundary into the temperate side at a decelerating rate, with . At any given time there is a wide range of x, above and below the transition temperature, where the temperature gradient is uniform. The sharp change of slope in the temperature curve lies at a temperature depression well below U _t at low resolution in temperature it is as if a freezing front were driving into an ice-water composite at 0¼. Since the corner occurs at Ζ = −0.22, the equation of motion of the freezing front is . The front moves from x = 0 to x = −14 m in 100 year and to x = −43 m in 1000 year. In the first hour it moves 14 mm.

  Two transition widths could be defined. The first, which we shall denote Δ₁Z, is defined as Δ₁ Z = U _t/U′(Z _t) and measures the range of Ζ in which the temperature changes by U _t, that is, the range in which the temperate-cold transition takes place. In our example

  

  which means 2.8 m after 100 year and 8.7 m after 1000 year. After lh this width is 3 mm. The second transition width, denoted by Δ₂Z, measures the width of the sharp break in slope seen in the U(Z) curve. This is suitably defined as , where is the slope at U = 2U(–∞); physically, it is the width in the temperate slab of the freezing front. In this example Δ₂Z = 0.003, which means 0.2 m after 100 year and 0.6 m after 1000 year. After lh this width is only 0.2 mm, which is less than the grain-size (about 2mm).

  The diffusivity of solute makes an insignificant difference for Ζ > 0. However, the freezing front is sufficiently sharp and travels sufficiently slowly for diffusivity of solute to have a marked effect on its width. The broken curves in Figure 7c are for D_c = 0. As we show analytically in Appendix C, the violent behaviour over a small range of Ζ is the result of a logarithmic dependence of X upon U.

Fig. 5. U: Ζ for the two-slab model with both slabs initially cold. The curve of C: Ζ would be almost indistinguishable from the curve draiun for U: Z. The asymptotic temperature depressions are indicated; these are the initial values for the two semi-infinite slabs.

Fig. 6. Curves of U, C and a/S^½ for the two-slab model with both slabs initially deeply temperate. Broken curves are for D_c = 0.

Fig. 7. Initially, the. lefthand slab is deeply temperate and the righthand slab is cold, (a), (b) and (c) are at successively higher magnifications. Broken curves in (c) are for D_c = 0.

(B) Solutions with C’(0) ≠ U’(0). Temperate: Temperate

This class of solutions differs from those just considered by being less smooth at Ζ = 0. Recalling that small Ζ corresponds to large t, we see that these solutions will tend to possess detail at large t resulting from diffusion processes slower than normal thermal diffusion, either anomalous thermal diffusion or diffusivity of solute.

A typical result is shown in Figure 8a. In this and the following graphs we again choose the impurity content at t = ∞ and the value of λ to be those of example (2) in section 2.5 (fine-grained glacier ice). Thus M (0) = 2.3 × 10^–8 and a = 6.0 × 10^–8. At Z = 0 (T = ∞) the relative values of U and C have been chosen so that nine-tenths of the temperature depression is due to concentration. This determines the temperature depression to be u = 2.3 × 10^–3 deg with r = 24 μm Thus U(0) = 1.41 × 10^–5, C (0) = 1.27 × 10^–5. This is well into the temperate range (U _t = 1.52 × 10^–4). The value of C′(0) has been deliberately chosen negative in contrast with the positive value of U′(0), and adjusted so that C (∞) C(–∞). Thus, the initial concentrations are about the same in the two slabs, but the temperatures and vein sizes are different. The curve for U is asymmetric and approaches U(–∞) faster than it approaches U(∞), because the effective thermal diffusivity is smaller on the low U (lefthand) side.

Fig. 8. Temperate: temperate transition with the initial values of C about the same in the two slabs. The detail near Ζ = 0 in (a) is resolved on a larger scale in (b). Broken curves are for D_c = 0.

There is no perceptible detail in U near to Ζ = 0, but this is in contrast to the curves for C and a/S ^½, which, although continuous, show rapid changes. These are better resolved in Figure 8b, which shows the same data at a larger scale. The interpretation is that, as already noted, there are two time-scales involved.

The corresponding curves for C and a/S ^½ when D_c = 0 are shown broken in Figure 8b. The curve for U is not shown because it is not perceptibly different from the full curve. The broken curves follow the full curves except near Ζ = 0, and at Ζ = 0 there are discontinuities in both C and S. This illustrates that, in the process with D_c ≠ 0, the mass of solute per unit length of vein in each slab remains nearly constant until the effect of D_c is felt at small Z.

As another example, Figure 9 has C′ (0) > U′(0). Here the initial vein size is larger on the colder side (Z > 0). This is achieved because initially the mass of solute per unit length of vein is 50 times larger on the cold side Table 1. The non-uniformity seen in the full curve of a/S ^½ has the form of a hump which is not centred on Ζ = 0. This indicates a travelling wave where the vein cross-section is constricted. (We recall that the graph represents a/S ^½ as a function of X at fixed T, and to visualize the development as Τ increases the whole graph is simply stretched uniformly, keeping the point at Ζ = 0 fixed, the stretching factor being T ^½) Since the maximum of the hump is at Ζ = −0.038, its equation of motion is . The wave moves from x = 0 to x = −2.3 m after 100 year and to x = −7.4m after 1000year. In the first hour it moves to x= −2.5 mm. It starts as very sharp and grows wider as it travels, but its peak amplitude remains constant. The existence of this wave is rather surprising; it simply emerges from the initial discontinuity between the slabs.

Fig. 9. Temperate: temperate transition with a travelling wave of pinched cross-section S on the warmer (lefthand) side. Broken curves are for D_c = 0.

In Figure 8, maxima and minima, not centred on Ζ = 0, appear in the curves for C; these represent travelling peaks and troughs of concentration, growing more diffuse as they travel but retaining their amplitudes. It is interesting that, if D_c = 0, the phenomenon of travelling peaks and troughs disappears.

7. Final Remarks

For simplicity, our computations were for uniform λ (length of vein per unit volume), although it would have been easy to allow two different values in the two slabs. Uniform grain-size is certainly never attained in glaciers and regions of different grain-size (different λ) may constitute one of the sharpest inhomogeneities that can persist.

To apply the results to glaciers, as was explained in section 1, it has to be remembered that the present model leaves out water flow along the veins (and along larger channels) and also the motion and deformation of the ice.

Therefore, conclusions about glaciers must be tentative. Whether a given transition thickness should be considered as sharp or broad depends, of course, on the context. The transition thickness of 200 m after 100 year for the cold cold example can be a significant fraction of the depth of a glacier and in that case the transition could not be modelled as sharp. On the other hand, some of the other transition thicknesses calculated, such as 2.8 m after 100 year for the deeply temperate: cold boundary, suggest that there are situations where a glacier model with an infinitely sharp boundary would be a useful approximation.