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 Reference ColbeckColbeck (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 Reference Brandt and WarrenBrandt and Warren (1993) to illustrate that the elevated sub-surface temperature in high-density Antarctic snow reported and modeled by Reference SchlatterSchlatter (1972) may be in error They refined the model of Reference SchlatterSchlatter (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 Reference SchlatterSchlatter (1972) and Reference ColbeckColbeck (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 (Reference Marshall and WarrenMarshall and Warren, 1987; Reference Brun, Martin, Simon, Gendre and ColeouBrun and others, 1989; Reference JordanJordan, 1991; Reference Greuell and KonzelmannGreuell 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. Reference Brun, Martin, Simon, Gendre and ColeouBrun 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 (Reference Brandt and WarrenBrandt 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 Reference JordanJordan (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 Reference Brandt and WarrenBrandt 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 (Reference Akitaya and ShimizuAkitaya and Shimizu, 1987) and wet snow (Reference IzumiIzumi, 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

Numerical Model

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

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 Reference PantankarPatankar (1980) and has been described in detail in Reference JordanJordan (1991). Thermal conductivity k_s is estimated from snow density as

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 Reference YenYen (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

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

Measurements of the downwelling longwave flux were unavailable and are estimated from a formula based on the equation of Reference IdSoIdso (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 (Reference Jordan, O’Brien and AlbertJordan and others, 1989). The turbulent fluxes are computed as (Reference Andreas and MurphyAndreas and Murphy, 1986):

and

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

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 (Reference Large and PondLarge 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 (Reference Grenfell and MaykutGrenfell and Maykut, 1977; Reference Brandt and WarrenBrandt and Warren, 1993) defined as

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 Reference Bohren and BarkstromBohren and Barkstrom (1974):

The empirically derived value of 0.003795 suggested by Reference AndersonAnderson (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

for the surface-control volume, and

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 Reference Bohren and BarkstromBohren 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 Reference Grenfell and PerovichGrenfell and Perovich (1984) for sparse clouds (case b) with spectral albedos estimated from the modified two-stream procedure of Reference Choudhury and ChangChoudhury 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 Reference WarrenWarren (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 Reference Choudhury and ChangChoudhury 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 Reference PerovichPerovich (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 Reference Grenfell and PerovichGrenfell and Perovich (1984) with spectral albedos computed from the formula of Reference Choudhury and ChangChoud-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 Reference Brandt and WarrenBrandt 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.