Sub-surface melting in a seasonal snow cover

Abstract

The ability of solar radiation to penetrate into a snow cover combined with the low thermal conductivity of snow can lead to a sub-surface temperature maximum. This elevated sub-surface temperature allows a layer of wet snow to form below the surface even on days when the air temperature remains sub-freezing. A high-resolution frequency-modulated continuous wave (FMCW) radar has been used to detect the onset of sub-surface melting in a seasonal snow cover. The experimental observation of sub-surface melting is shown to be in good agreement with the predictions of a one-dimensional mass- and energy-balance model. The effects of varying snow characteristics and solar extinction parameters on the sub-surface melt characteristics are investigated using model simulations.

Introdution

The study of the physical processes that occur near the surface of a seasonal snow cover is complicated by diurnal variations in the surface-energy balance, resulting primarily from swings in air temperature and solar radiation. This was illustrated by Colbeck (1989), who modeled the temperature profile in a seasonal snow cover by assuming a sinusoidally varying surface temperature and exponential decay of solar radiation within the snow cover. His calculations showed that a sub-surface temperature maximum can occur, if the rate of solar heating beneath the snow surface is greater than the rate of cooling by conduction. A more refined energy-transport calculation at the snow surface was conducted by Brandt and Warren (1993) to illustrate that the elevated sub-surface temperature in high-density Antarctic snow reported and modeled by Schlatter (1972) may be in error They refined the model of Schlatter (1972) by using a narrow spectral resolution for solar radiation to take into account the wavelength-dependent extinction in snow. The narrow-band model indicated that the broadband approximation used by both Schlatter (1972) and Colbeck (1989) overestimates the magnitude and depth of the sub-surface maximum and that, for clean snow, significant sub-surface temperature elevations are limited to low-density snow covers. Other researchers (Marshall and Warren, 1987; Brun and others, 1989; Jordan, 1991; Greuell and Konzelmann, 1994) have used two- or three-band models that offer significant improvement over the one-band variety without the computational inefficiency of fine spectral resolution. Brun and others (1989) observed the occurrence of temperature maxima a few centimeters below the snow surface on cold, sunny days and successfully replicated measured temperature profiles using a three-band model divided at wavelengths of 0.8 and 0.15 μm.

The modeled results of the elevated sub-surface temperature effects have been difficult to verify experimentally. This is primarily due to the difficulties involved in measuring the temperature profile of a snow cover using temperature sensors embedded in snow Radiative heating of these sensors on clear days can lead to incorrect temperature measurements (Brandt and Warren, 1993) More conclusive evidence of elevated sub-surface temperatures may be the onset of sub-surface melting, which can be easily detected using a high-resolution frequency-modulated continuous wave (FMCW) radar. FMCW radars have the characteristics of detecting layers within a snow cover, if sufficient electromagnetic contrast exists between the layers, arising from discontinuities in dielectric constant, snow grain-size and/or surface roughness. The large dielectric contrast between ice and water at the radar Geqnencies makes it possible to detect subsurface melting with an FMCW radar.

We have recently observed that sub-surface melting is not an uncommon event in a seasonal snow cover. In this paper, a case study is presented of a sub-surface melting event observed with an FMCW radar. A one-dimensional mass- and energy-balance model for snow (SNTHERM), developed by Jordan (1991), is used to predict the temperature and liquid-water content profiles of the snow cover. The solar absorption routine in SNTHERM uses a simple two-band resolution of the solar spectrum (0.4–1.12μm and 1.12–2.4μm), in which the longer-wavelength component is “dumped” at the surface and the shorter-wavelength fraction is absorbed exponentially. The experimental observations of sub-surface melt are shown to be in good agreement with the modeled results. These encouraging results gave us confidence to run additional simulations in order to test the effects of varying snow characteristics and solar-extinction parameters on the sub-surface melt characteristics. The simulations indicated that the occurrence of sub-surface melt is strongly dependent on snow density and generally concurred with the observation of Brandt and Warren (1993) that sub-surface temperature elevations are minimal in high-density snow.

The ability to detect and model the sub-surface melt in a seasonal snow cover can contribute to a better understanding of the physical processes that occur near the snow surface. Practical applications of these subsurface effects may be found in snow remote-sensing and avalanche investigations. Microwave interaction with a snow cover is greatly influenced by the presence of liquid water; therefore the conditions that can lead to an increase or decrease in snow wetness are of considerable interest to the radar remote-sensing community. Prediction of avalanche releases may be improved by investigating the processes that contribute to the formation of mechanically weak snow layers. Avalanche releases have been attributed to layers of depth hoar (Akitaya and Shimizu, 1987) and wet snow (Izumi, 1989) which were buried by subsequent snowfalls. It is possible that subsurface melting and subsequent freezing of this melted layer can also affect the mechanical strength of the snow below the surface.

Experiment

a. Environmental conditions

A case study of the FMCW radar observation of subsurface melting that began around 1230 h on 1 March 1993 in Hanover, New Hampshire, is presented. The subsurface melting was observed on a calm, clear day when the air temperature was below freezing. Figure 1 presents some of the relevant meteorological conditions during the onset of melting. The air temperature remained sub-freezing for most of the day except between 1400 and 1800 h (the maximum temperature of 2.8°C was observed around 1600 h). During thia period, the hourly averaged wind speed ranged between 1 and 3 ms⁻¹ The incoming solar radiation measured with a pyranometer indicated that the sky was relatively free of clouds during the onset of sub-surface melting. The fractional cloud cover reported at an airport 6 miles [9.6 km] to the south confirmed the measurements of the pyranometer.

Fig. 1. Hourly summaries of the meteorological conditions during the onset of sub-surface melting.

The snow-cover properties were obtained from snow pits located approximately 5 m from the radar footprint. The snow cover consisted of several layers whose total depth was approximately 30–40 cm (due to the uneven contour of the radar plot, there was some spatial variability). Profiles of the snow temperature, density and grain-size observed at 0900 h are reported in Table 1. The thermal conductivity and effective thermal conductivity are computed from algorithms within the computer simulation to be discussed later. The snow temperature was obtained by inserting at varying depths a single thermocouple probe that was attached to a digital thermometer. The temperature of the entire snow cover was sub-freezing and had a negative gradient due to overnight radiative cooling of the surface.

Table 1. Initial snow-cover characteristics and control-volume grid used in computer simulation for 1 March 1993

The top 4 cm of the snow cover fell on 22–23 February and by 1 March had densified to 130 kg m⁻³. Below the top layer was the remains of 20 cm of snow that fell on 16–17 February, which had compacted to approximately 15 cm and had a density of 170 kg m⁻³. The snow grain-size of the top two layers was approximately 0.5 mm. The bottom third of the snow cover consisted of a 2 cm crust on top of approximately 10–20 cm of snow. This layer was created by the snowfall of 12–13 February, when approximately 20 cm of snow fell over an existing shallow snow cover. During the latter part of this storm, the snow was wet and formed the crust when it subsequently froze. The density of the bottom layer was approximately 250 kg m⁻³.

b. FMCW radar measurements

The FMCW radar system that was used to detect subsurface melting is illustrated in Figure 2. An HP 8350B sweep oscillator and an HP 83554A millimeter-wave source module generates signals whose frequency varies linearly with time from 26.5 to 40 GHz with a sweep time of 80 ms. A directional coupler divides the frequency-modulated signal into two paths, a reference signal is brought directly into the mixer (directional coupler and crystal detector) and a target signal is routed to the snow cover via a transmitting antenna. The reflected signal from the snow cover is fed into the mixer input via a receiving antenna. An HP 3562A signal analyzer is used to obtain the fast Fourier transform (FFT) of the mixer output. The FFT yields a power spectrum where the frequency is proportional to the return travel time (function of the target distance and the effective dielectric constant of the propagation path) of the source signal and whose magnitude is proportional to the target reflectivity. An example of a power spectrum, illustrating the output of an FMCW radar, is illustrated in Figure 3. The first peak in the power spectrum at 6.012 kHz represents the reflection from the snow surface. The second peak at 6.248 kHz represents the reflection from the crust approximately 19 cm below the surface. For the FMCW radar system used in this study, the relations between the distance within the snow, ΔΖ(cm),and the frequency difference in the power spectrum, Δf(Hz), is expressed as

(1)

where t_swp is the sweep rate (80 ms), bw is the band width of the FMCW sweep (13.5 GHz), c is the velocity of light and ϵ_eff is the effective dielectric constant of snow. For a record length of 80 ms, the signal analyzer has a frequency resolution of 12.5 Hz, which corresponds to a resolution of approximately 1 cm in low-density snow.

Fig. 2. Schematic of the FMCW radar. DC and CD are directional coupler and crystal detector, respectively.

Fig. 3. FMCW radar output. The frequency of the power spectrum is proportional to the electromagnetic distance and the magnitude is proportional to the target reflectivity.

For this study, the FMCW radar was mounted on a gantry approximately 6 m above the ground (the distance between the base of the antenna and the ground was 5.6 m) and positioned so that the transmitted signal was directed toward the snow cover at a 4° incident angle. The radar back-scatter from the same spot on the snow cover was recorded once every 30 min. A sequence of FMCW radar returns illustrating the formation of a water layer beneath the snow surface is shown in Figure 4. The return at 1210 h is very similar to that observed at 0910 h (cf. Fig. 3), indicating that very little change in snow property occurred during the 4 h period between these two traces. The formation of a melt layer beneath the snow surface began to emerge around 1240 h until a distinct layer was detected around 1310 h. The frequency difference between these two peaks was 0.025 kHz. For snow, whose density is 130 kg m⁻³ (effective dielectric constant of approximately 1.28), this frequency difference indicates that the top of the melt layer occurred approximately 2.0 cm below the surface (the thickness of the melt layer cannot be inferred from the radar measurement) Application of water-sensitive dye to the snow below the surface confirmed that a layer of wet snow due to sub-surface melting was present.

Fig. 4. Sequence of FMCW radar response illustrating the formation of a melt layer below the snow surface

Further analysis of the radar data indicated that the maximum melt appeared to have occurred around 1400 h. The increase in liquid water resulted in higher absorption loss, so that the magnitude of the reflection from the crust was at its minimum. After 1440 h, the magnitude of the reflection from the crust layer gradually increased, which indicated that the melted layer had begun to freeze. The reflection from the crust gradually increased and then stabilized after 1740 h, which indicated that liquid water was no longer present. Comparison of the FMCW radar returns from 0910 and 1740 h shows that, due to the sub-surface melt and freeze cycle, the FMCW radar was able to detect a refrozen layer that was not present before the melt.

The ability to detect sub-surface melting with an FMCW radar has been demonstrated. The FMCW radar provides information about the time and location of the melt layer. To investigate the physical processes responsible for the sub-surface melt, the radar data are analyzed using a one-dimensional heat- and mass-balance model.

Numerical Model

Conservation of energy for a one-dimensional snow soil system is described by the partial differential equation

(2)

where c_s is the effective heat capacity of the wet snow, ρ_s is the snow density, T is temperature (Kelvin), M is the rate of phase change (or melt) from ice to liquid, L_m is the latent heat of fusion, t is time, k_e is an effective thermal cunductivity, F_solar is the net solar radiation (positive upwards) and z is the vertical position measured upwards from the snow-soil interface. Converted heat resulting from water flow is not included in Equation 2, since the liquid-water content in this application does not reach the minimum level necessary for flow. Sensible- and latent-heat effects are combined in an apparent heat capacity, which incorporates a freezing curve for snow artificially selected so that it approaches a step function at 0°C but is not too steep as to pose difficulties with the numerical solution. The numerical scheme is based on the control-volume methodology of Patankar (1980) and has been described in detail in Jordan (1991). Thermal conductivity k_s is estimated from snow density as

(3)

where k_a, and k_i, are the thermal conductivities of air and ice, respectively, and the coefficients have been selected so that Equation 3 fits the data of Yen (1962) and extrapolates to the cunductivity of ice. Latent-heat effects due to the diffusion and sublimation of water vapor are incorporated in the effective thermal conductivity as

(4)

where L_v is the latent heat of sublimation, D_e is the effective diffusion coefficient of water vapor through snow and ρ_v is the water-vapor density.

The surface-energy flux F_top, (positive upwards) combines the net solar flux (F_solar, the net longwave-radiation flux (F_ir) and the turbulent-energy fluxes of sensible (F_sen) and latent heat (F_lat), giving

(5)

Measurements of the downwelling longwave flux were unavailable and are estimated from a formula based on the equation of Idso (1981). As is customary, the emitted component of the upwelling longwave flux is computed from the Stefan-Boltzmann equation, in which the snow emissivity is taken as 0.97 (Jordan and others, 1989). The turbulent fluxes are computed as (Andreas and Murphy, 1986):

(6)

and

(7)

where ρ_a is the air density, c_a is the specific heat of air at constant pressure, C_H is the bulk-transfer coefficient for sensible heat, Τ_a is the air temperature, Tⁿ is the temperature of the surface-control volume, n, C_E is the bulk-transfer coefficient for latent heat, w is the wind speed, ρ_v ⁿ is the water-vapor density of the surface-control volume and ρ_v,a is the water-vapor density in air, computed from the relative humidity. The bulk-transfer coefficient for sensible heat is computed from the roughness length Ζ₀, giving for neutral stability

(8)

where Z′ is the observation height above the snow interface. A value of 0.007 m is used for z₀, The bulk-transfer coefficient for latent heat at neutral stability C_EN is taken as 0.7C_HN. For unstable atmospheric conditions, the standard stability adjustment is made to C_EN and C_HN (Large and Pond, 1982) but, for stable conditions, the correction is omitted, since it was found to degrade the performance of the model significantly.

The magnitude of the solar-absorption term depends primarily on the intensity, angle and spectral composition of the incident radiation, and on snow albedo, optical depth and grain diameter. As a first approximation, incident solar energy is assumed to be diffuse and isotropic. The heating rate is computed using a net bulk-extinction coefficient β_ext (Grenfell and Maykut, 1977; Brandt and Warren, 1993) defined as

(9)

where β_λ, is the spectral flux-extinction coefficient, F_λ (0) is the net spectral radiation at the snow surface and h is the snow depth. The spectral integration is simplified into two broad bands corresponding roughly to the wavelength ranges of 0.4–1.12 μm and 1.12–2.4 μm. Radiation within the lower-energy band is assumed to be totally within the top-control volume. That within the higher-energy band is assumed to decay according to Beer’s law, using a net asymptotic bulk-extinction coefficient β_∞. computed from the formula of Bohren and Barkstrom (1974):

(10)

The empirically derived value of 0.003795 suggested by Anderson (1976) is used for the adjustable parameter C. Combining the contributions from the two energy bands, the discretized energy gain due to solar heating within the snow cover is

(11)

for the surface-control volume, and

(12)

for the interior control volumes j, where F^j+½ is the net solar radiation at the upper boundary of control volume j and f is the fraction of net solar radiation in the higher-energy band. Since heat fluxes are defined as positive in the upward direction, a negative value for the change in solar radiation ΔF_solar corresponds to absorption of energy by the snow cover. Although the fraction f will vary with cloud cover and snow characteristics, a constant value 0.46 is used here, which corresponds to the best fit between predicted and estimated surface temperatures for the standard 10 d test period used in validating the model.

Model Simulations

Temperature and liquid-water profiles in the snow cover were simulated for 1100–1900 h on 1 March 1993, which contained the period of melt. Values of snow density and grain-size used in the simulation are provided by the profiles shown in Table 1. Within the two upper layers, a control volume thickness of 0.5 cm was used. The total snow depth of 34 cm was sufficiently deep that the appeared semi-infinite to solar radiation. The simulation was begun at 1800 h on the previous day to allow the temperature profile to come to equilibrium prior to the period of interest. The observed temperatures at 0900 h were used to verify the simulation and compared favorably. Values for the thermal and effective thermal conductivity at 1100 h are shown in Table 1 and were computed from Equations (3) and (4).The thermal conductivity of the upper layer is relatively low, typical of low-density snow, and latent-heat changes due to the diffusion and sublimation of water vapor account for about 50% of the effective thermal conductivity. A summary of the surface-energy components for the simulation period a shown in Table 2, computed from meteorological observations recorded at 15 min intervals and fractional cloud-cover observations from a nearby airport. The spectrally integrated albedo for 1 March (computed as the ratio of upwelling to downwelling measurements of solar radiation) is plotted in Figure 5 and averaged 0.81 for the day. A rise in albedo in the morning generally corresponded to a period of increased cloud cover, when the spectral composition of the incident radiation was shifted towards the shorter wavelengths.

Fig. 5. Spectrally integrated albedo of snow (ratio of upwelling to downwelling measurements of solar radiation from Eppley pyranometers). 1 March 1993, Hanover, New Hampshire, U.S.A.

Table 2. Summary of the surface-energy components for the simulation period

Hourly predictions of volumetric liquid-water content are shown in Figure 6. The simulation generated output on a 15 min basis but for illustrative purposes only the half-hourly values are plotted. Sub-surface melt is first predicted at 1215 h at a depth of about 2.75 cm and reaches a maximum of 2.06% at 1630 h at a depth of 4.25 cm. Although the maximum coincidentally occurred at the top of the underlying denser snow layer, a later run on homogeneous snow (described below) showed a maximum at the same depth The timing of melt agrees closely with that of 1240 h observed with the radar. The top of the simulated melt layer occurred at about 2.75 cm depth, somewhat below the observed depth of 2.0 cm. At around 1545 h, the calculated total volume of liquid water in the snow cover reaches a maximum, about 1 h later than the observed time of 1440 h (when the radar return from the ice layer was at a minimum). Temperature profiles on the hall-hour are shown in Figure 7. A maximum surface temperature of −l.5°C is reached at 1430 h. As can be seen in Figure 8, the heat gained by the top-control volume from the external sources of radiation and turbulent transfer from the atmosphere is negative. Therefore, any warming of the snow surface must be through conduction of heat from the underlying region of elevated temperature, which in turn is warmed by solar absorption. After sunset, the snow surface cools rapidly, dropping to −9°C by 1830 h. Later in the afternoon, the melt layer a calculated to freeze, both from the top down and the bottom up, and is completely refrozen by 1845 h. By comparison, radar observations showed the water to have refrozen by 1740 h.

Fig. 6. Hourly predictions of volumetric liquid-water content. Simulation is initialized with the snow characteristics in Table 1.

Fig. 7. Simulated temperature profiles on an hourly basis.

Fig. 8. External energy flux affecting top-control volume (net longwave, radiation sensible- and latent-heat fluxes, solar radiation absorbed by top-control volume).

Several alternative simulations were run to test the effect on melt depth of varying snow characteristics and extinction parameters. The following parameters were varied: the high-energy fraction f, the constant C in the Bohren and Barkstrom (1974) bulk-extinction equation, snow density, grain-size and roughness length. The depth, timing and magnitude of maximum melt for these runs are summarized in Table 3. In order to test clearly the effect of density, the snow cover was taken as homogeneous for the sensitivity studies.

Table 3. Depth, timing and magnitude of maximum melt for variations in f, C, snow density, grain-size and roughness length

Trial 1 is the base-line simulation and is the same as the case described above, except that the snow characteristics are taken uniformly as those of the top layer. Accordingly, the density is 130 kg m ⁻³ and the grain-size is 0.5 mm. The timing and depth of maximum melt are the same as for the layered simulation but the magnitude is less because less energy is absorbed by the uniformly lighter snow (recall that previously the maximum has occurred at the top of the older snow layer). In trials 2 and 3, the fraction of radiation “dumped” in the top-control volume is varied. An approximate fraction of net radiation in the 1.12–2.4 μm band was computed by convolving the incident spectral-irradiance data of Grenfell and Perovich (1984) for sparse clouds (case b) with spectral albedos estimated from the modified two-stream procedure of Choudhury and Chang (1979). In the integration procedure, a spectral resolution of 0.005 μm was used and the indices of refraction for ice were taken from Warren (1984). The resulting fraction of 0.53 agreed closely with the value of 0.54 used in the baseline case. Increasing the energy dumped in the top-control volume to 0.75 resulted in no melt and decreasing it to 0.25 markedly increased the amount of sub-surface liquid water.

The net bulk-extinction coefficient is varied in trials 4 and 5. For the base-line case, trial 4 and trial 5, B_∞ computes as 22.06, 11.63 and 58.13 m⁻¹, respectively. Increasing the extinction coefficient (trial 5) resulted in a shallower melt depth and an earlier occurrence of the peak. The depth of the melt layer observed with the radar was best replicated using a value of about 40 m⁻¹ for Β_∞, A further trial was conducted in which the two-band routine in SNTHERM was replaced with a simplified four-band model, in order to test whether a narrower spectral resolution would improve the simulation. The spectrum was divided into bands of 0.4–0.93, 0.93–1.12, 1.12–1.39 and 1.39–2.4 μm and effective bulk-extinction coefficients were computed from Equation 9, using the depth (h−z) at which half the radiation within the band was absorbed. Spectral albedos and spectral bulk-extinction coefficients used in the integration were estimated from the asymptotic formulas of Choudhury and Chang (1979). The resulting estimates of β_∞ for the four bands (high to low energy) were 7.50, 33.46, 77.29 and 516.64 m⁻¹. The depth of the upper edge of the melt zone for the four-band solution was unchanged in comparison with the base-line case but the maximum melt occurred at the slightly shallower depth of 3.75 cm. While the simulation was marginally improved, refining the spectral resolution did not appear to account for all the residual error. A plausible explanation for the discrepancy between theory and observation is that the snow at the Hanover site was slightly contaminated, leading to somewhat higher extinction coefficients than expected for pure snow. Because the natural snow cover was inhomogeneous, it is also possible that a multi-layered model, such as that described by Perovich (1989), is necessary accurately to model sub-surface melt.

Varying the grain-size by itself (trial 6) had the same effect as changing the parameter C, since both parameters affect the magnitude of B_∞ . A computed all-wave albedo of 0.83 was close to the observed value of 0.81, giving us confidence that the estimated grain-size of 0.5 mm was reasonable. If the site were contaminated, however, it is possible that the actual grain-size was smaller. The all-wave albedo was computed by integrating the spectral irradiance data of Grenfell and Perovich (1984) with spectral albedos computed from the formula of Choud-hury and Chang (1979), again using a 0.005 μm spectral resolution.

When the snow density was increased to 250 kg m ⁻³ (trials 7 and 8), no melt occurred. Increasing the density increased both the extinction coefficient and the thermal conductivity but the thermal conductivity appeared to influence the occurrence of sub-surface melt more strongly. Because of the density-dependence in Equation 3, the effective thermal conductivity was approximately 2.5 times that of the base-line case and the region of elevated temperature was more rapidly cooled by conduction. The sub-surface temperature for trial 7 was still elevated (about 0,7°C above that of the surface) but never reached the melting point.

The depth of the melt zone was moderately sensitive to changes in surface roughness, as illustrated by trials 9 and 10. The surface roughness governs the magnitude of turbulent-energy exchange between the snow cover and the overlying air. Since the combined turbulent fluxes are negative until 1430 h, larger and smaller values of z₀ will result in greater and lesser cooling οf the surface, respectively. Correspondingly, the melt zone shifted upward and downward For trials 9 and 10. The radar observations are best matched with the shorter roughness length of 1 mm, both in the location and earlier timing of the melt. In addition to uncertainties in the computation of the turbulent fluxes, the estimated downwelling longwave flux may have been in error due to the approximate cloud-cover values used in the simulation.

Conclusions

A high-resolution FMCW radar has been used to detect sub-surface melting in a seasonal snow cover. A one-dimensional mass- and energy-balance model using a simple two-band resolution for the solar spectrum is shown to be effective in reproducing the experimental results. Additional model simulations were run to investigate the effects of varying snow characteristsics and solar-extinction parameters on the sub-surface melt characteristics. In agreement with Brandt and Warren (1993), sub-surface melt in the model simulations was limited to low-density snow. It is reasonable to conclude that on calm, clear days when the air temperature is near freezing, sub-surface melting is likely to occur in low-density snow. Since these conditions are not uncommon in a seasonal snow cover, the importance of sub-surface melt and its effect on the near-surface snow properties and processes should not be minimized.