1. Introduction
One approach to the problem of a polythermal glacier (Reference Fowler and LarsonFowler and Larson, 1978; Reference HutterHutter, 1982; Reference FowlerFowler, 1984; Reference Hutter, Blatter and FunkHutter 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 (Reference MaderMader, 1990, Reference Maderin press a, Reference Maderb). 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 (Reference LliboutryLliboutry, 1971; Reference HarrisonHarrison, 1972; Reference Nye and FrankNye and Frank, 1973; Reference Raymond and HarrisonRaymond and Harrison, 1975; Reference Harrison and RaymondHarrison 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.
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 103 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 (Reference WorsterWorster, 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 (Reference NyeNye, 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 (Reference Raymond and HarrisonRaymond 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 (Reference MaderMader, 1990, Reference Maderin 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 (Reference BearBear, 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 Reference RobinRobin (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
where vi, is the specific volume of ice, sw 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 (Reference SchwerdtfegerSchwerdtfeger, 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 2 = αr v 2.
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 Reference Walford, Roberts and HillWalford 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 (Reference NyeNye, 1991a) has revealed a systematic error in the method arising from diffraction, and new measurements of the diffraction pattern (Reference Walford and NyeWalford and Nye, 1991) have given ϕ = 25.0 ± 1.0° for a supposedly typical grain boundary. At the same time, there is evidence (Reference MorrisMorris, 1972; Reference MaderMader, 1990, Reference Maderin 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 Reference Nye and MaeNye and Mae (1972) to α = 0.1007.
Equation (1) may now be written, using r, as
where
With νi = 1.09m3Mg−1, γ = 34mJm−2, sw – s i = 1.22 kJ kg−1 deg−1, we find A = 5.44 × 10−9m 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−2. 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
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
where Dth = 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 Dc 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 2λ. So the flux per unit area of continuum is –D cπr 2 λ∂c/∂x. Since the mass of solute per unit volume of continuum is πr 2λc, this leads to the conservation equation
We generalize this for veins distributed is otropically to
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 cos2θd λ. Since the factor cos2θ 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 Reference MorrisMader (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 Dc = 0.
Then Equation (6) gives
, and, if we denote by m the mass of solute per unit length of vein (m = constant), we have
Equation (2) is now
where the constant Β = Κm/π.
Since u is now a function only of r, and vice versa, the second term in Equation (4) involving ∂(r 2)/∂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,
or
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
where
in agreement with Reference HarrisonHarrison (1972, p. 26). ut is a transition temperature depression. When u falls to ut, the specific heat capacity becomes twice its normal value. This simplified formula is valid when u ≫ A 2/B or u ≫ πA 2/Km.
The behaviour of the effective diffusivity Deff 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)Dc, 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 Dc rather than Deff that then controls the approach to equilibrium, as we shall show by numerical solutions in section 6. We shall call such ice “deeply temperate”.
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 Deff < Dc/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 (Reference MaderMader, in press a,Reference Maderb), was observed to have r = 6 µm at a temperature depression u = 0.10 deg. In example (2), from Reference Raymond and HarrisonRaymond 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 3 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−4 deg for example (2), which we neglect in comparison with the measured u. Then B = ur 2 = 3.6 × 10–12 deg m2 for (1), and 1.2 × 10−12deg m2 for (2). Taking λ = 5 × 104 m−2 for (1) and 1 × 106 m−2 for (2) and L/σ = 160 deg, gives ut = 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 2/B = 1.6 × 10–5 deg for example (1) and 5.0 × 10−5 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 Deff falls to (1/3)DC (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
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
where just two dimensionless constants remain, namely,
With D c = 1.3 × 10–9 m2 s–1 and D th = 1.2 × 10–6 m2 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
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
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).
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
Equation (7) is simply
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 Dc → 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 Dc is some 103 times smaller than the diffusivity for heat in ice, one might be inclined to put Dc = 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.
Α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−4), 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 Dth dominates when |Z| ≈ 1, while Dc is effective when |Z| ≈ 10–3. The numerical schemes adopted are discussed in Appendix A.
To illustrate the role of Dc, computations were also made (Appendix B) with the same asymptotic values as before, but setting Dc = 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.
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−8, a = 6.0 × 10−8. 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 Reference HarrisonHarrison (1972), according to whether its temperature depression is below or above the transition temperature depression ut 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), ut 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),
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 Dc 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>>Ut, 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 Dc = 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−3). The lefthand slab is deeply temperate with initial temperature depression U(—∞) = 1.03 × 10−6. 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−4, U′ (0) = 3.40 × 10−3).
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 = Ut
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 Δ1Z, is defined as Δ1 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 Δ2Z, 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 Δ2Z = 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 Dc = 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.
the final temperature depression is the mean of the two initial temperature depressions.
, 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.
, is 4.6 times less than U(0). So the common boundary in this example is cold. 
, or, in terms of physical variables,
. 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.
, where 
is the slope at U = 2U(–∞); physically, it is the width in the temperate slab of the freezing front. In this example Δ2Z = 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).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.
Curves of U, C and a/S½ for the two-slab model with both slabs initially deeply temperate. Broken curves are for Dc = 0.
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 Dc = 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.
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 Dc = 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 Dc = 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 Dc ≠ 0, the mass of solute per unit length of vein in each slab remains nearly constant until the effect of Dc 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.
Temperate: temperate transition with a travelling wave of pinched cross-section S on the warmer (lefthand) side. Broken curves are for Dc = 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 Dc = 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.
Acknowledgements
I am grateful to Dr M. Ε. R. Walford and Dr H. Mader for many helpful discussions on the behaviour of veins in ice, and for kindly commenting on the paper. The experimental work by Dr H. Mader, undertaken to test the hypothesis referred to in section 2.3 about the behaviour of the impurities as the veins freeze, was supported by a grant from the U.K. Natural Environment Research Council.
The accuracy of references in the. text and in this list is the responsibility of the author, to whom queries should be addressed.
APPENDIX A Numerical Methods
To integrate Equations (17), (18) and (19) outwards from Ζ = 0, two different schemes were used. The first was suitable for small |Z|, where a finer interval was often needed. However, for larger Ζ this scheme was sometimes unstable and then the technique used was to switch over to the second scheme.
The first scheme was as follows. Differentiating Equation (17) and using the new variable E(Z) = U′(Z), the Equations (17), (18) and (19) become
The starting values at Ζ = 0 for U, U′, C, C′ yield corresponding starting values for E,E′,C,C′,S,S′. Then Equations (A1), (A2) and (A3) are used to integrate forward for S(Z), E(Z) and C(Z).
For larger Z, Equations (17), (18) and (19) were used directly. Equation (19) was solved for S′:
and Equation (17) is
Knowing U,C,S at two successive points labelled i – 1 and i, Ui+i is obtained from Equation (A5), Ci+1 is computed from Equation (A6) in the form
and Si+1 is obtained from Equation (A4).
APPENDIX B Solutions withDc = 0
We consider here the solution of Equations (17)–(19) in the approximation Dc = 0, or equivalcntly d = 0. Equation (19) then integrates, with Equation (21), to give
Thus M, the mass of solute per unit vein length, which is initially uniform in the two slabs, now remains uniform and constant and so independent of Z. However, in general, it will have different values in the two slabs and we call these M+ and M–. Equation (17) now becomes
where M now stands for M + or M–.
To integrate the remaining Equation (18) we need S′. The relation
combined with dU/dS from Equation (B2), converts Equation (18) to
Equation (B2) may be solved as a quadratic for S in terms of U, and hence Equations (B3) has the form
This is solved by outward integration from Ζ = 0 with trial values (7(0) and U’(0), the appropriate values of M+ and M– having been computed in advance from Equation (B1) using prescribed values of C and S at Ζ = ±∞. This trial integration will result in asymptotic values (U(∞), U(-∞) which are different from those required. Therefore U(0) and U′(0) are now adjusted, by trial and error, until the correct asymptotes are obtained. We now have curves for Dc = 0, and for Dc set at its correct value, which have the same asymptotes; that is, they refer to the same initial temperatures and impurity contents in the slabs. Note that the final temperature depression U(0) is different according to whether solute diffusion is allowed or not.
APPENDIX C Logarithmic Solution
Figure 7c shows extremely sharp corners in the broken curves for Dc = 0, the transitions occurring over a Ζ interval of order 10−3. We examine this behaviour analytically.
The graphs show that, in the region concerned, the major contribution to the temperature depression is made by the solute concentration. So, as an approximation, place a = 0 in Equation (23) to give
and then Equation (B3) reduces to
a differential equation for U. Since we are well into the temperate region, for which the transition temperature depression is M½, by Equation (22), the first term in parentheses in small. Furthermore, as a local approximation wc may put Ζ as constant, Ζ =-Z0 say (Z0 > 0). Then Equation (C1) becomes
Introducing a rescaled temperature depression
reduces this to
which is readily integrated, first to give
and then again to give
where p is a constant. As Z → –∞, V → 1/p from above. The other constant simply represents an arbitrary translation in Ζ and, because the solution is only valid near Ζ =-Z0 , we place it equal to Z0. The asymptote is low down on the V-axis for large values of p, and this is the case wc want. As V increases by 1/p Ζ changes by order 1/p 2. Thus, the lower is the asymptote the sharper is the corner. In our example Ρ = 40.
