1. Introduction
The dynamic role of the liquid-water content in temperate ice remains one of the unsolved problems in glaciology despite the fact that it has been shown to affect flow at the melting point (Reference LliboutryLliboutry, 1976; Reference Lliboutry and DuvalLliboutry and Duval, 1985) and that first theoretical models (Reference HutterHutter, 1982; Reference FowlerFowler, 1984; Reference Hutter, Blatter and FunkHutter and others, 1988) have pointed out its significance in the global boundary conditions. In fact, very little is known about water content in temperate glaciers. Typical values found by Reference Vallon, Petit and FabreVallon and others (1976) in an ice core of an alpine glacier are between 0 and 4%.
Two methods have been proposed for measuring liquid water in glaciers, (i) The dielectric constant of the mixture ice-water depends on the relative presence of the two components and may serve as an indicator of water content. But this method, used for snow, is not sensitive enough and the dielectric constant may vary with other parameters such as grain boundaries or point defects, (ii) The calorimetric method consists of measuring the amount of energy needed to freeze the water in a given volume of temperate ice. It has the advantage of being very sensitive, and the results are less influenced by the other parameters. However, the use of an adiabatic calorimeter requires the extraction of an ice core. Stresses in the probe are thus changed, and this modifies the water content (Reference RaymondRaymond, 1976).
The need for a calorimetric method for in-situ measurements (i.e. in the bulk of the glacier, without extracting a core) led us to study the problem of the propagation of a “cold wave ” in an ice—water mixture. The wave front consists of a surface of phase change that separates cold ice from the ice-water mixture. Its speed depends on the imposed boundary conditions on the specimen and the water content. By determining the migration of this well-defined boundary, we may thus infer the water content in the specimen (Fig. 1).
Fig. 1. Temperature of the ice and position of the interface during the experiment.
The method of determining w by this procedure is an inverse problem, because its distribution must be inferred essentially from information gained along the boundary of the domain. Such problems tend to be difficult and sensitive to measurement and numerical error. Here, however, two simplifications may be justifiably invoked which make the problem tractable and, as our results show, reasonably stable to numerical error. First, length scales that are considered are a few centimetres (less than 0.25 m) over which water contents in the vein system and the inclusions may be regarded as constant. Secondly, the region outside the phase-change surface may be considered to be a “bath” with uniform conditions. This bath will only affect the solution of the freezing problem at the phase-change surface. This implies that the water content should only vary very slowly over length scales of the above-mentioned 0.25 m. Field observations lead us to be confident about these assumptions.
Approximate analytical solutions of the solidification problem exist for one-dimensional geometries (see Reference VujanovicVujanovic, 1989). We chose to solve the problem numerically as it seemed most suitable for the various boundary conditions we consider.
2. The Initial Value Problem
2.1. Hypotheses
Consider a specimen of temperate ice at 0°C. Assume that a metal mass has been inserted at its center and that the contact between the metal and the temperate ice is perfect. We may consider two types of Gedankenexperiments:
-
Assume that at time t = 0 the metal mass is suddenly set at a prescribed temperature T0 below the freezing point and the system is subsequently left free to evolve.
-
Assume that at time t = 0 the metal mass is subject to a temperature jump and that its new temperature (below freezing) is kept constant while the surrounding ice is free to adjust its thermal state to the new conditions.
In both cases, the cold will propagate from the metal into the ice and freeze the interstitial water. The cold front that forms the surface of phase change between the wholly cold ice and the ice-water mixture can be identified experimentally as the location where the temperature gradient suffers a jump and thus a kink in the temperature profile arises. Observations have shown these to be clearly identifiable.
We impose the following idealized initial conditions:
The temperate ice is an homogeneous mixture of water and ice uniformly at the melting temperature = 0°C.
The specimen is thought of as infinitely large so as to eliminate the influence of finite boundaries. The cold source (metal) is an infinite sheet, an infinitely long wire, or a spherical ball, i.e. linear, axial, or spherical symmetry prevails.
The metal is a perfect conductor and the contact between the ice and the metal is perfect. (This condition will be relaxed later on.)
2.2. Equations
Because of the spatial symmetry, only positions with 0 < r < R(t) need be considered. The metal mass ranges from r = 0 to r = a, the ice from r — a to r = m, whereas the cold ice occupies the domain a ≥ r ≥ R(t), R(I) being the position of the phase-change surface. In this latter domain, the heat-conduction equation is
where ρt is the ice density (= 918 kg m−3), Ci is the heat capacity of ice (= 2.12 χ 103Jkg−1K−1), ki is the thermal conductivity of ice (= 2.2 W K−1 m−1), Θ is the temperature, and ∆ is the Laplacian. Throughout this text, the dot represents the time derivative.
Boundary conditions have to be satisfied at the interface “cold ice-temperate ice ” r = R(t) and at r = a. At the former, the heat flow from the phase-change surface is balanced by the latent heat released by freezing, viz.
at r R(t), t > 0, and
Here, L: is the latent heat of fusion, w is the moisture content just ahead of the freezing front, and ∂/∂n denotes the spatial derivative perpendicular to the phase-change surface.
At r = a, heat flow and temperature must be continuous. However, for very small values of a and a perfect conduction, we may regard the metal as a body with uniform temperature and heat capacity VcρcCc,. Here, Vc is the volume of the metal (copper) piece. The time rate of change of the heat stored in the copper, Vcρcccθ must then be balanced by the heat flow through the contact surface with the ice,
at r = a, t > 0. This condition applies to Gedanken-experimenl (1). When the temperature is held constant, we have
at r = a, t > 0. Finally, we have the following initial conditions:
To simplify the equations, we use the problem symmetries and then introduce dimensionless variables. If the cold source is a plate, a cylinder or a sphere, we will use Cartesian, cylindrical, or spherical coordinates. Equation (1) will be respectively
Equation (2) can be written as
and Equation (4) becomes
where λ = 1, 2, 3 for Cartesian, cylindrical, and spherical geometry, for which equation numbers are (4′), (4″), and (4‴), respectively. Next, we introduce the dimensionless variables r = Xox, I = t0t′, θ = u0u, R = X0S, where
for experiment (1) when the initial heat content is prescribed, and
when, as in experiment (2), the temperature of the metal is held fixed. The problem can then be summarized by the equations given in Table I.
2.3. Numerical solution
We describe in this paragraph a numerical solution for the equation with cylindrical symmetry at fixed energy. The difficulty induced by the moving boundary can be avoided by use of the boundary position instead of the time as second independent variable, as shown byGarofalo (unpublished), which introduces an artificially fixed domain of integration and is possible because the position of the phase-change surface is a monotonically increasing function of time. We shall then define
TABLE 1. Non-dimensional equations
where ϕ is the speed of the phase-change surface. As s(t) is a strictly growing function, the time corresponding to a given position is calculated by
With the new variables u′, s, and ϕ, the equations for the cylinder become (upon omitting the primes)
for a ′ < x < s and a′ < s < ∞, subject to the boundary conditions
at x =a′, and
along the line x = s.
The numerical scheme for its solution is described in Appendix A.
3. Variable Coefficients And Insulating Layer
In this section, we introduce the thermal variability of the thermodynamical quantities. We also examine the effect of an insulating layer between the ice and the cold source. Only the case of a constant temperature source with cylindrical symmetry will be developed.
3.1. Variable coefficients
Allowing for variations of the coefficients, Equation (I) becomes
Equations (2), (3), and (4) remain the same as before. The dimensionless variables are identical to those of section 2.1 with the mere addition of
where ρi ci and ki are reference values at 0¼C, and primes denote dimensionless quantities. Replacing the independent variable t by the position of the moving boundary, we obtain the initial boundary-value problem (primes are again omitted)
for a′ < x < s,
for x = s and a < s < ∞, and
for x = s = a′, of which a solution scheme is constructed in Appendix B.
Initial values are the same as before. In the calculation we use for к and pc the values given by Reference Hobbsobbs (1974, p. 360–61)
where
is the thermal diffusivity.
3.2. Insulating layer
We will now model the influence of an insulating air layer between the cold source and the ice. The heat transfer in the layer is supposed to be stationary. In this case, the heat flow through a layer of thickness ε is, for a cylinder of unit length, given by
η being the thermal conductivity of the layer. The boundary condition at x = a now becomes (Reference Carslaw and JaegerCarslaw and Jaeger, 1959, P- 19)
The introduction of the dimensionless variables defined by
gives a system formed by Equations (13), (14) and
at x = a′; see Appendix B for the numerical peculiarities.
4. Numerical Results
4.1. Test of the numerical solution
An analytical solution for the freezing of water exists in the case of a constant temperature T0 at x = 0 (classical Stefen problem slightly modified for an ice-water mixture with a water content iv). The cold-ice temperature T is given by
in which all variables are dimensional and where ξ is given by
D is the thermal diffusivity of the ice, с is the heat capacity, and L is the latent heat of melting. The position of the interface in time is described by the equation
This equation allows us to test the algorithm. The results presented in Figure 2 show that for large values of the dimensionless initial velocities, C1 = d.s/d/(0), the phase-boundary position is accurately predicted, the error being of the order of 0.1% when x < 0.01 [m]. Other tests with different values of r, λ, hx, and hs (for definitions see Appendix A) have demonstrated that the most suitable values are r = λ = 0.5 and hx = hs = 0.01. However, a problem arises for very short time.
Fig. 2. Error between analytical and numerical solutions for constant coefficients and mirror symmetry: G1 = (ds/dt)(0).
In a similar way, the equation having temperature-dependent coefficients has been tested by blocking ils temperature dependence and comparing the results with runs for the constant-coefficient algorithm. The results show again a satisfactory large time behaviour with an important divergence for short time. Nevertheless, this solution can be trusted for t > 1 [s].
Fig. 3. Time for a displacement of the boundary of 10 cm as functions of T0/w (in 0 °C). Solution at constant coefficients for the cylinder with fixed temperature and fixed energy, the plate and the sphere with fixed energy.
Fig. 4. Phase-boundary velocity plotted against its position. Onset of the oscillations for the experiment at fixed temperature, as shown in the inset. T0/w in °C.
4.2. Constant coefficients
As it appears from our scalings, we can for the case of constant coefficients combine the two parameters T0 and w, and present the results as functions of T0/w. The time required for a phase-boundary movement of 10cm is given in Kigure 3 for the different geometries of Gedankenexperiment (I) and for the cylindrical symmetry of Gedankenexperiment (2). Alternatively, Figure 4 shows for the experiment with a fixed initial heat content the phase-boundary velocity as a function of its position. The inset displays the continuation of the graph with T0 = −50°C when oscillating “instabilities” arise. These oscillations are typical of Crank-Nicholson schemes and cannot be avoided, but they falsify somewhat the evaluation of the travel time of the phase boundary. The onset of such oscillations is expected whenever the available energy is very near the energy needed for the boundary displacement. It can be seen in Figure 3 that the oscillation begins at different values of the parameter T0/w for the plate, the cylinder, or the sphere, owing to the different amounts of energy that are needed per unit displacement of the phase boundary. The difficulty is due to the use of the phase-boundary position instead of time as an independent variable. This method breaks down when the phase boundary approaches a static equilibrium position. Under these conditions, another
Fig. 5. Time for a displacement of 10 cm plotted against T0/w (in C) for several values of w. Solution for variable coefficients for the cylinder.
computational technique would be more suitable (see. for example, Reference CrankCrank, 1984).
In the experiment at fixed boundary temperature, the computation yields convergent results at all values of T0/w (Fig. 4).
4.3. Variable coefficients
With the variable-coefficients scheme, results depend on TQ/w and w ; many more iterations are required for the convergence test to be successful and their numbers increase with increasing values of |T0/w|. The maximum number of iterations is as high as 11 at the beginning of the computation, compared with 4 for the constant-coefficient equation. Moreover, the larger w is, the more will the dimensionless coefficients k and ρc in Equation (13) vary and thus destabilize the scheme, as has been observed by Garofalo (unpublished). But for reasonable values of T0 and W, the solution is satisfactory.
For comparison with Figure 3, we display in Figures 5 and 6 the time for a boundary displacement of 10cm as a function of T0/w and for different values of w. Comparing Figure 5 with Figure 3, we infer that, as expected, the displacement for larger values of iv is more rapid if we take into account the increase of thermal conductivity with a decrease of temperature. Also of interest is the enormous impact of an insulating layer of air as small as 0.1 [mm], which tends to hide the water-content effect (Fig. 7).
Fig. 6. As Figure 5 with an insulating layer 0f air of 1 mm.
Fig. 7. Position of the boundary as a function of time for a fixed temperature of -5 ° C and various water contents.Solution with variable coefficients.
Figure 7 shows the phase-boundary position versus time for an initial temperature of −5°C and different water contents, in the absence of an insulating layer. It can be seen that the measurement of the boundary position, as an indicator of the water content, should yield a good sensitivity, and therefore be suitable as a measuring technique.
5. Conclusion
An analytical solution for solidification has been described for various types of boundary conditions. The problem with constant coefficients gives satisfactory results. The variation of the coefficients introduces an instability at high values of the water content of temperate ice, w. Nevertheless, the scheme works for small values of w. In the view of a practical application of this problem to the measurement of the water content, we infer that this method should provide reasonable sensitivity. However, a small insulating layer between the cold source and the ice causes an important delay in the movement of the phase boundary and tends to hide the influence of the water content.
Appendix A
We outline a numerical solution scheme for Equations (9)-(12). Supposing that a solution u(x,s) exists, we can discretize the domain as shown in Figure 8 and write
Fig. 8. Fig. 8. Discretization of the domain.
Using finite differences, Equation (9) can now be replaced by the approximate equation
where r and λ are two weights used to stabilize the calculation and hx and hs are the step sizes. This particular scheme has been used to avoid problems arising from the initial conditions. We need for this equation the value of Uj +1,j which lies outside the computational domain. This value'is approximated by using a Taylor series expansion for the function u, as will be explained below. After rearrangment, we obtain
with
and
Moreover, condition (10) at x = a′ can be written as
Equations (A1) and (A3) form a tridiagonal system which can be solved for the unknowns ui j+1 i=1 …, j, if the values of u¡ j and q>y are known. The system is first solved for Γ j = (h X/h S)φ j , φ j+1 is then updated (we shall explain below how this is done) and the system solved again with Γj given by Equation (A2), until the following convergence test, proposed by Garofalo (unpublished), is verified:
where N is the number of iterations.
To complete the solution scheme, we still need an equation to update φ j + 1,
the values at the boundary uj j and uj+1, j,
the initial values u 1,1, u 1,2, t 1, t 2, φ1, φ2,
an equation to calculate the time tj+1
Using a Taylor series expansion for u on the boundary X = s, we obtain
From Equations (9) and (11) we deduce
Given that on the line x = s
and that dx/ds =1, we infer ∂x/∂s = -∂u/∂s and thus obtain, in view of Equations (9) and (11)
dx2x
on x = s. Thus, with condition (12)
In much the same way, we derive
and Equation (12) implies
for j > 1. At time t = O when x = s = а′, a combination of Equations (10) and (II) yields
Furthermore, as is evident from Table I, and from Equation (A3) with j = 1 we have
With Equation (A2), Equation (A4) for j = 1 yields φ2. By integrating Equation (8) with Simpson’s r ule, we finally have
where φ3/2 is given by linear interpolation between ψ1 and φ2
We now have all the necessary ingredients for a numerical solution of the stated initial boundary-value problem in cylindrical coordinates when the initial heat content in the metal specimen is given. We first compute u1 j with Equation (A7a), then obtain If1 from Equation (A6), U1 2 from Equation (A7b), and φ2 from Equation (A4). With these quantities being determined, Γ1 follows from Equation (A2), and u¡ 2 (i = 1,2, …) can be computed from Equations (Al), (A3)', and (A5). The ensuing steps then only involve Equations (Al), (A2), (A3), and (A5) with the iteration suggested on Γ and Equation (A8) to update the time. A solution for other symmetries is easily derived in the same way. In the case of the problem at fixed temperature, only the boundary condition (10) is changed, and Equation (A3) is replaced by
Appendix B
We construct here a numerical solution scheme for Equations (13)-(15). To solve Equation (13) numerically, we adopt the finite-difference scheme of Garofalo (unpublished, p. 17) and set
and
Re-arranging and using Equation (14), Equations (13) and (15) imply the following tridiagonal system
The boundary conditions at x = s are obtained by writing Equation (13) in the form
and by employing the same reasoning as in section 2.2, explicitly
From these relations we deduce in much the same way as before
When the equations for the insulating layer are solved (section 3.2), the new tridiagonal system is composed of Equations (Bi) and
Equations (B3)-(B6) and (B7) also apply here, but the initial conditions are
Equations (B3)-(B6) and (B7) also apply here, but the initial conditions are
from Equation (16), and
from Equations (14) and (16).
The system can then be solved with the same algorithm as before.
