Introduction
The growth of sintered bonds between ice spheres has been studied (Reference Nakaya and MatsumotoNakaya and Matsumoto, 1954, Reference JensenJensen, unpublished; Reference Hosier and HallgrenHosier and Hallgren, 1960; Reference KingeryKingery, 1960; Reference KuroiwaKuroiwa, 1961; Reference Hobbs and MasonHobbs and Mason. 1964). Reference Hobbs and MasonHobbs and Mason (1964), concluded from the time-dependence of the growth of the bonds and from the measured absolute value of the diffusion constant that bond growth is caused primarily by the transport of water molecules through the vapour phase. However, the strength of bonds in the initial stage of sintering (after short times of contact between the grains) cannot be explained by this transport mechanism.
The strength of snow is largely dependent on intergranular structure and on the strength of the bunds between the grains. Deformation of snow often causes fracturing of some of the bonds and the generation of new contact points. Deformation changes the intergranular structure (Reference KryKry, 1975[a] and Reference Kry[b]; Reference Lawrence and BradleySt. Lawrence and Bradley. [1975]; Reference GublerGubler, 1978; Reference LawrenceSt. Lawrence. 1980). The relationships between stress, strength and strain, and strain-rate depend strongly on the strength development in the new contact points.
In order to improve our knowledge of the initial stage of sintering, we designed an experiment to measure the strength of bonds between ice grains after contact times of I to 1000 s as a function of contact pressure, temperature, and stress rate. The experimental results are explained by the flow of water to the bond neck in the liquid-like surface layer of ice (Reference WeylWeyl, 1951; Reference FletcherFletcher. 1968; Kvlividze and others, 1973; Reference Golceki and JaccardGolecki and Jaccard, 1977; Reference GilpinGilpin, 1979), the deformation of the spheres at the contact point, freezing of the liquid-like layer in the contact volume, and local pressure melting.
Strength increase of snow caused by sintering for short time periods (1–1000 s) at new contact points is the basis for deterministic models for locking mechanisms in avalanches, blowing snow, snow compaction, and snow milling.
Apparatus
The basic experimental arrangement is shown in Figure 1. The tips of two ice cones are slowly brought into contact and are compressed using forces which range from 0 to 0.1 N. One of the ice cones (diameter of the base: 20 mm, point angle 45°, tip curvature radius 0.5 mm) is mounted on a piezoelectric force transducer (sensitivity < 10−3 N). The second, moveable cone is fixed to a steel–aluminium axle which is suspended by thin threads. The forces (in compression and in tension) are applied to this rod by electromagnetic coils. The system is placed in a sealed container. The container bottom is partially covered with granular snow and can be heated with a temperature-controlled heater to maintain the desired relative humidity of between 90 and 100% in the container. The experimental area with the two cones is in addition protected by a small “Plexiglas” box. Temperatures are measured with an absolute accuracy of ±0.1 °C at different points around the cones and inside the ice cones near their tips (Fig. 1). The polycrystalline ice cones were slowly frozen in brass moulds from their base to the tip under a temperature gradient using outgassed, distilled water. The electronic arrangement achieves an automatic soft contact of the two cone tips and can apply a large variety of stress functions in tension and compression to the contact between the cones.
Basic experimental setup.
Measurements
Recordings of bond strength were made at temperatures between −0.2 °C and −44 °C. Measurements with two different pairs of cones were made at each temperature. After the cones were carefully brought into contact, a uniform force pulse of duration t and amplitude 10−2 N was applied in compression. The compression pulse was immediately followed by a tensile force applied at a rate of 2 × 10−1 N s−1 until the bond fractured. The cones were kept apart for 60 s between consecutive experiments to allow the cone tips to recover. The duration t of the compression pulse was varied according to the following pattern:
alternately forward and backward. Each run consisted of at least three patterns of the compression pulse lengths and included one change of the cone pair. Alternation of the start of the runs between the compression pulses of the longest or the shortest duration is important to test for any possible hysteresis in the recorded tensile fracture forces.
Results
The recordings of the fracture forces type A (Table 1) were plotted as a function of the compression pulse duration t for each run and fitted to a power law:
α and β are temperature- and geometry-dependent variables. The power law was chosen because all of the currently proposed growth laws for sintered bonds between ice grains (Reference Hobbs and MasonHobbs, 1974) can be expressed as power laws of the time of contact of the ice spheres. For each run three pairs of parameters α and β, their mean errors, correlation coefficients, and t-values (Student’s t-test) were calculated, α↑,β↑ for increasing compression pulse length t, α↓, β↓ for decreasing t and the mean values ᾱ, β for the complete run.
The following estimates were used for the determination of the relative error of α and β:
where σ F is the mean error of the smoothed fracture force F T . The determination of the compression pulse duration t is assumed to be error free. t min = 2 s and t max = 320 s are the lower and higher limits of t. β(1 ± r β ) corresponds approximately to the lowest and largest slope angles β in a log plot of the fracture force versus compression pulse duration t (Fig. 2). The mean correlation coefficient for all fits is 0.9, the extreme values are 0.5 and 0.99. The corresponding confidence levels are always higher than 95%. Results for β and α are plotted in Figures 3 and 4 as a function of the negative temperature δ[|°C|]. The mean relative errors are r α = 0.22 – 0.28 (range for α↑, α↓, ᾱ ), r β = 0.17, the corresponding most probable errors arer α = 0.13 –0.2,r β = 0.06 − 0.11. The error distributions are not symmetric but are skewed to larger relative errors caused by enlarged errors at high temperatures.
Fracture force versus compression pulse length t at a temperature of −10 °C. The straight lines are the power law fits for increasing and decreasing t and the whole run.
β↑ as a function of the negative temperature δ[|°C|] fitted with a seventh order polynomial. The corresponding curves for β↓ and β are also included.
α↑ as a function of δ fitted with a logarithmic function of δ for δ < 1 deg and a third degree polvnorninal δ ≧ 1 deg.
The relationship between the fracture force and the compression force is roughly linear (measurements type B, table I). The gradient is about 1, independent of temperature (see page 471). All measurements show a residual fracture force for zero contact pressure.
The fracture force does not significantly depend on stress rate (measurements type C, Table I).
Discussion of the results
The hysteresis of the fracture force at increasing and decreasing compression pulse durations is caused by surfaces being incompletely recovered. No distortion of the cone tip surface could be seen under an optical microscope a few seconds after the separation of the cones. The hysteresis of the mean fracture force is equal to or less than the mean error of a single fit. The mean relative hysteresis of the exponent β is 0.15. Above −10°C the higher exponents are related to measurements with increasing compression pulse durations t, the lower exponents to recordings at decreasing t. Below −10 °C no significant hysteresis has been observed.
The parameter β has its highest value in the temperature range −1°C to −5°C, falls to a slightly lower constant value at lower temperature but above −1 °C decreases sharply reaching zero at the melting point. The fracture-force proportionality parameter α increases sharply above −2°C. Near the melting point, the fracture force reaches a high asymptotic value after short contact times (1–30 s). The increase of the fracture force with contact time is highest at temperatures of about −1 °C to −5 °C. Below −10 °C the magnitude of the tensile fracture force decreases only slightly with decreasing temperature; the growth rate of the bonds remains almost constant.
Theory
Assumptions of the model
The main assumption in the model is the existence of a liquid-like layer on ice surfaces. This liquid-like layer has been proposed by several authors (Reference FaradayFaraday, 1859; Reference WeylWeyl, 1951; Reference KingeryKingery. 1960; Reference FletcherFletcher, 1968; Reference Jellinek and HorneJellinek, 1972; Kvlividze and others, 1973; Reference Golceki and JaccardGolecki and Jaccard, 1977; Reference Maeno and NishimuraMaeno and Nishimura, 1978; Reference GilpinGilpin, 1979). Fletcher is the only worker who has tried to develop a quantitative theory of the surface structure of ice; both his attempts have been shown to be incomplete. His theories predict surface thicknesses of the order of magnitude of some 10−9 m with sharp increases as temperature approaches the melting point. Reference Jellinek and HorneJellinck (1972). using mechanical testing, estimated the surface thickness to be 9 × 10−8 m at −1.8 °C. Reference Golceki and JaccardGolecki (1977) measured the thickness of the disordered surface as a function of temperature using proton channelling. Measuring the thickness of the transition layer as a function of temperature shows a logarithmic dependence on the absolute value of the negative temperature in |°C| (≈log (const. δ[°C])) similar to that stated in Fletcher’s original theory (Reference FletcherFletcher, 1962) but with a 20-fold larger slope. Golecki’s measurement is also in fair agreement with Jellinek’s mechanical estimate. The nuclear magnetic resonance measurements of Kvlividze and others (1973) show a significantly increased mobility of a small percentage of the molecules in dispersed ice. The mobility was found to be comparable to the mobility of molecules in liquid water. Kvlividze has shown that the high mobile molecules are concentrated at the ice–gas interface.
With these results as background, I propose the following model: If two spherical ice particles are in contact, the liquid-like layers start to freeze in the contact area. Surface energy together with the large difference in the curvatures between the neck and surface points on the spheres outside the contact area give rise to a pressure gradient along the transition layer. This pressure gradient caused by the surface tension creates in the transition layer a viscous flow toward the bond neck. It is further assumed that the liquid-like layer retains its equilibrium thickness in the source region independent of the flow. The flow distance L is the mean distance at which liquid contributes to the bond. L is estimated in different ways: (a) to be a constant value; (b) as the distance between points of minimum and maximum pressure; (c) as a function of a relaxation time describing the phase change at the solid–liquid interface. If the particles are compressed, the spheres deform plastically until the stresses in the contact region decrease below the yield stress. In the initial stage of contact under compression, pressure melting contributes significantly to the bond volume at high temperatures.
Bond geometry
The idealized parameters are defined in Figure 5. The curvature of the neck surface is approximated by a circle of radius r K.
where
-
r b = bond radius in the neck midplane
-
R = radius of curvature of the ice particles in the contact region
-
d = penetration depth of the particles because of plastic deformation.
Bond geometry.
The radius of the bond at the neck boundary r* is
The volume V′ of the sphere cut is
The total bond volume V is given by the following expressions
The curvature k b in the mid-plane cross-section of the bond is
The curvature of the particle in the contact region is
Driving pressure gradient
The pressure difference Δp over the flow distance L is approximated by
γ LV is the surface energy of the transition layer.
Under the assumption of a laminar Newtonian flow with a constant pressure gradient parallel to the surface layer and constant flow height determined by the equilibrium conditions for the disordered surface layer, the rate of volume increase dV/dt of the bond volume Vis (contribution of both particle surfaces considered)
where l is the equilibrium thickness of the liquid-like layer and η L is the viscosity of the surface layer.
Flow distance L
The estimation of the flow distance is critical. The dependence of the bond radius on the contact time and particle diameter is partially determined by the approximation of the flow distance. The simplest assumption is a constant flow distance L, independent of geometry and material flux (case a):
As we will see later, only the now distance multiplied by the viscosity of the surface layer can be estimated from the experiment.
In a second attempt, the flow distance is chosen to be twice the distance over which the driving pressure gradient is established by the surface energy and the bond geometry (case b). For this case the flow distance depends on the actual bond geometry and the grain diameter:
A third possibility is to assume a constant relaxation time at the ice–liquid interface
where i is a typical intracrystal molecular distances≈5 × 10−10 m and τ is the relaxation time which controls the phase change between the ice structure and the liquid structures ≈10−3 s (Reference Bilgram and GüttingerBilgram and Giittinger, 1978).
The flow rate
where A is the source area and
The flow length is now a function of the flow rate:
Creep under contact stress using the indenter -hardness concept
The yield stress σ y for creep of polycrystalline ice is given by Hobbs (1974)
with E = 2.6 × 10−19 J, n = 3.9 for T≥ 263 K, E = 1.2 × 10−19 J, n = 4.4 for T< 263 K = Boltzmann constant: 1.38 × 10−23, J deg−1, and t is the time of compression; the constant c has the dimensions [s1/n ].
σ 0 is determined by experiment and is of the order of 106 N m−2 for t = 1 s, T = 270 K. For ihe times t after the application of the compression force, the reciprocal penetration of the cone tips is calculated as a function of the corresponding yield stress:
r′ is the effective radius of the stress-conducting part of the bond and F is the compression-force pulse amplitude.
The bond radius in the absence of material flux r c is approximately
The compression force F c acting on the cross-sectional area πr 2 c is deduced from the total compression force F with the assumption of equal yield stress for the cone ice and the bond ice:
From Equations (15) and (16) the penetration d is
The corresponding contribution to the total bond volume is (Equation 5)
Pressure melting
The depression of the melting point ∆T M as a function of pressure is
where T M is the melting point of ice under normal conditions = 273.15 K, v L, v s are the specific volumes of the liquid and the solid phases. L LS is the specific latent heat of fusion, and C≈ 7.4 × 10−8 K m2 N−1.
The equilibrium penetration depth d is reached if ΔT M is equal to the negative particle temperature δ in °C.
The corresponding contribution to the total bond volume is similar to Equation (18).
Heat flow
The heat flow from the growing bond volume (refreezing of the liquid) to the different sources of the liquid (pressure melting, surface flow) limits the possible equilibrium growth rate of the bond volume. The supercooled liquid in the bond volume refreezes at a speed proportional to its supercooling. For a rough estimate, we assume that the liquid refreezes from the sphere’s surfaces into the bond volume. The maximum possible refreezing distance is equal to h (Fig. 5).
The radius of the refrozen bond is approximately
With a linear refreezing velocity v G = gδ (g is the factor of proportionality and equals 10−4 m s−1 deg−1, Reference Bilgram and GüttingerBilgram and Güttinger, 1978) we get an approximate value for the corresponding heat flow Φ by multiplying an approximate bond surface area by the refreezing velocity and the specific latent heat L LS(= 3.35 × 108 J m−3):
for V ⪢ 0 and r s ⪢
The heat flow between source and sink is roughly
where b is the distance between source and sink, k is the heat conductivity (= 2 J m−1 s−1 deg−1). Δ is the temperature difference between source and sink, and A is the cross-sectional area for heat flow. A is of the order of magnitude of the bond cross-sectional area πr 2 b, b is approximately equal to L (flow distance). Equations (22) and (23) yield for the temperature difference δ between source and sink
For the initial states of sintering with r b/R ⪡ 1, L ⪡ R, and (r b/R)L < 10−4 m, the temperature difference Δ between source and sink is always small compared with δ (temperature below the melting point in °C) for δ > 0.1 deg and contact times > 2 s.
Thickness of the liquid-like layer
Reference FletcherFletcher (1962) and Reference Golceki and JaccardGolecki and Jaccard (1977) propose a logarithmic dependence of the transition layer thickness l on supercooling:
The revised model of Reference FletcherFletcher (1968) and the Golecki and Jaccard measurements indicate an even sharper increase of the layer thickness above −1 °C. Our model calculations agree with the measured data only when the log (thickness–temperature) relationship is modified in the same direction (Fig. 6).
Thickness of the transition layer as a function of temperature δ.
This dependence is in rough qualitative agreement with Fletcher’s revised model and is in fair quantitative agreement with the measurements of Reference Jellinek and HorneJellinek (1972) and Reference Golceki and JaccardGolecki and Jaccard (1977).
Approximate solutions
Basically, Equation (11), which describes the rate of increase of the bond volume, has to be solved for the bond radius r b. Assuming a temperature- and time-independent strength of the bond ice, the cross-sectional area πr 2 b would be proportional to the measured tensile fracture force.
Unfortunately, the time-dependent variation of the penetration depthd (Equation 15) and the complicated expressions for the bond volume V (Equation 6) and the flow distance L (Equation 13), as well as the additional effects of pressure melting and of the transition layer impede a closed form solution of Equation (11).
However, the following approximate solutions can be found:
With the assumption d ⪡ R we get the following estimate for r b
Using this approximation, the differential Equation (12) can be solved for
where
Rewriting liquation (12) for L constant gives:
For L ≈ πr k:
For
where
-
r K ≈ r b 2/2R (Equation 3),
-
r b ≈ (RV)1/4 (Equation 27),
-
r* ≈ r b + r K (Equation 4).
For these very limiting assumptions we get for the bond radii
and
The corresponding tensile fracture forces are
where σ 1 is the tensile strength of the bond ice for the given stress distribution as a function of contact time and temperature.
The time and temperature dependence of σ 1 is not known, therefore, we assume σ 1 to be independent of the age of the ice as well as of temperature. For the tensile fracture force the power of the contact time t varies for the different models for the flow distance L between β = 0.29 and 0.43. The power of the grain radius R varies from 0.71 to 0.86. The numerical solution will show that the effective values are lowered by the plastic deformation process. For t = l s, γ Lv = 0.1 Jm−2, η L = 2.5 × 10−3 kgm−1 s−1, R = 5 × 10−4m, δ = 1 deg, l = 8 × 10−8m, L(const.) = 3μ (Equation 30a). i = 5 × 10−7 τ (Equation 30c), the bond radius varies from 55 μm to 75, μm. With an ice bond strength of 106 N m−2 at a stress rate of about 106 N m−2 s−1 the tensile fracture force is in the range of 10−2 N (α in Equation 1). The comparatively low bond ice strength σ 1 (about 1/2 of the tensile strength measured by Reference ButkovichButkovich (1954)) may be caused by the high surface to cross-section ratio of thin bonds, the low age of the bond ice (some seconds old), and the complicated inhomogeneous stress state in the bond cross-section.
At low temperatures the main process is plastic deformation of the cone tips. From Equations (16 and 17) we have
The decrease of the power of t with decreasing temperature from 0.3–0.4 (viscous flux regime) to 0.23 (plastic deformation regime) is the main reason for the measured decrease of β (Equation 1, Fig. 3). Equation (32), according to my measurements, indicates that for the plastic deformation regime the fracture tensile force is proportional to the compressive force.
In reality, the two regimes are always mixed, plastic deformation being more important in the initial state of sintering and at low temperatures. If there is no compressional force, bond growth is completely described by the viscous flow regime for short contact times. For contact times longer than 102–103 s the main transport mechanism is vapour transport with a bond area growth proportional to t 0.4 (Reference Hobbs and MasonHobbs and Mason, 1964).
Numerical solution
The flow diagram for the simulation program is given in Figure 7. The program calculates r b(t) incrementally for given temperatures.
Flow diagram for the numerical solution.
Results
Plots of r b(t) at different temperatures (F = 10−2 N), σ 1πr 2 b(T) for t = 1 s, F = 10−2 N; σ 1πr 2 b (F) for t = 20 s, T = 263 K; β(T) for F = 10−2 N, 1 s < t < 300 s are presented in Figures 8–11 and compared with the experimental data.
Calculated bond radius rb(t) for temperature between −0.1 °C and −40 °C and a compression pulse amplitude of 10−2 N.A: total bond radius. B: radius of refrozen bond.
Tensile fracture force FT = σ 1πr2 b as a function of temperature δ(2) for t = 1 s, F = 10−2 N. Included are the measured data for ᾱ.
Calculated fracture force FT as a function of the compression amplitude F for contact times t = 20 s, T = 263 K. The band shows the range of the measured data (Table I, type B).
Calculated and measured values of β(δ) for compression forces F = 10−2 N and compression pulse lengths 1 s < t < 1 000 s. Measured data β (Table I, type A); calculated data.
2: σ y = 106 N m−2, 3: σ y = 1.8 × 106 N m−2, 4: σ y = 4 × 106 N m−2, 5: Flux regime Equation (30c).
Equation (13c) was used for the model calculation presented here (constant relaxation time at the interface between transition layer and bulk ice). The parameters are given below:
Correspondence between the experimental data and the theory is at its worst at temperatures above −5 °C. The maximum theoretical value for the bond growth rate is below −1 °C. Its value of c. 0.35 (simulation) to 0.4 (flux regime approximation; Equation 30c) is in fair agreement with the experiment. The theoretical value for the growth rate drops off faster with decreasing temperature than does the experimental value. This transition from the flux regime to the plastic deformation regime depends mainly on the temperature dependence of the yield stress (n, E in Equation 14). An increased yield stress at temperatures above −5 °C would drastically increase the correspondence between theory and experiment (Fig. 11). Additional, more precise measurements are necessary in this temperature range.
The sharp decrease of β(T) (Equation 1) at temperatures approaching the melting point is caused by the decreasing pressure gradient in the transition layer for bond radii r b approaching its saturation value. Because of the large transition layer thickness, bond growth is very fast for temperatures above −1 °C. However, at temperatures near the melting point, the bonds remain partly unfrozen for contact times up to several seconds (Fig. 8).
The tensile strength of the bond ice is assumed to be independent of temperature. Reference ButkovichButkovich (1954) measured a 25% linear increase in tensile strength between 0 °C and −40 °C. In addition, an ageing effect for increasing t possibly increases strength for higher contact times. Both effects, which are not included in the model, tend to increase correspondence between simulation and measurements.
Conclusion
The model describes the development of sintered bonds between ice grains for short contact times and compressive forces in fair agreement with the measurements. Improved measurements of the tensile fracture force F T and of the bond radius and direct recordings of the deformation d are necessary to refine the model.
Acknowledgements
I should like to thank R. A. Schmidt and K. C. Agrawal for their careful review of the manuscript.
