Introduction
There have been many studies in which heat flux in the surface snow layer is estimated, for example Dalrymple and others (1966) and Weiler and Schwerdtfeger (1977). But, there seem to be no investigations on the errors involved. Since the estimated values of snow heat flux are used for heat budget considerations on the snow surrace, error estimates and methods of correction are essential.
The Japanese Antarctic Research Expedition (JARE) carried out a micrometeorologica! observation at Mizuho Station (70°42'S, 44°20'E, 2230 m) as a part of Polar Experiment (POLEX-South) in Global Atmospheric Research Programme (CARP) (Mae and others 1981). Data compiled by Wada and others (1981) include hourly values of snow temperature, measured with platinum resistance thermometers at 0 (surface), and depths of 0.1, 0.3, 0.5, 1, 3 and 5 m, from 23 February to 31 December 1979. A main object of the project was to determine the heat budget on the Antarctic snow surface. This vast amount of snow temperature data made it possible to investigate the calculation of snow heat flux. In this research, the errors are obtained as a function of interval of estimating the flux, and correction methods are proposed for the major errors in short term estimation. The former part of this research is described in detail by Kikuchi (1983), thus the main part of this paper is devoted to the correction.
Heat Flux Calculation Formula
Snow heat flux at the surface can be calculated from numerical differentiation and integration of the well known divergence theorem ρc(?T/?t) = - ?S/?z, where T is the snow temperature, S the heat flux p the snow density, c the specific heat of ice, t the time, and z the depth (positive downward). The lower boundary is taken at a depth of 10 m where S is assumed to be zero and T is assumed to be constant at the mean temperature at 5 m depth.
The actual formula for calculation of the heat flux at the surface S0 becomes:
where zi is the depth of snow temperature observation, and the quantities with suffix i represent those at that depth. Time is represented by the second suffix t, and 2Δt is the time interval of flux estimation. The temperature is averaged over 2Δt before the calculation of Equation I in order to reduce errors described later. The specific heat of ice c(T) is 2.0 x 10-3 J kg-1 K-1 at -26°C, and is a slightly varying function of T. Thesnow density ρ is approximated as:
from data by Narita and Maeno (1978) on the snow/ice density at Mizuho Station. The discontinuity at 8 m in Equation 2 is within the observed scatter in the data.
Errors in Heat Flux Calculations
Four factors are considered as causes of errors in calculating heat flux with Equation 1:
-
a) Errors and/or resolution in the temperature measurement (0.1 K), especially in determining the temperature derivative in Equation 1.
-
b) Errors in the numerical integration, which assumes a straight-line profile between two levels of temperature observation.
-
c) Errors in the bulk snow density as determined by Equation 2 (30% error at the surface and decreases with depth to 0.5 % at 10 m).
-
d) Errors associated with the lower boundary condition (a rough estimate gives an error of 2 x 10-1W m-2).
Figure 1 shows the errors arising from factors a) c) and d) as a function of the estimation interval 2Δt. The total of the three effects above are also shown in the figure. The error from factor a) decreases rapidly as the estimation interval is lengthened because of the averaging process. But it occupies a major part of the total error, which is as large as 90% when the daily values are calculated. Errors by c) and d) are not so significant but limit the improvement at longer interval to about 10%.
Fig. 1. Errors in snow heat flux estimation as a function of interval of estimation. Temperature is averaged over a period which is equal to the interval.
The error from factor b) becomes large if the assumption of straight line profile is invalid. A measure of validity is obtained by comparing the layer thickness d to a characteristic diffusion length defined by 
where K is the thermal diffusivity and ω is the radian frequency of temperature variation. Table 1 lists 1 values for the layers of present consideration. The values of K are assumed on the basis of Yamada and Hasemi (1974) and Dalrymple and others (1966).
If d is larger than 1, Equation I yields considerable error. The values listed in Table 1 suggest that only slow variations can be estimated from this measurement. The total error in daily heat flux shown in Figure 1 is further increased by the factor b).
table. 1. Characteristic Length For Temperature Variations
One of the other factors not included here is the direct heating of the sensors by shortwave radiation. The values of surface temperature were compared with those of air temperature which had been measured under shielding. The correlation was fairly good and there was no evidence of excessive heating in the surface temperature. Therefore, we did not take the radiation error into account, but more precise experiment would be required for testing the radiation effect.
Correction For Temperature Error
As shown in Figure 1, the error caused by the temperature resolution (a) is significant in short term heat flux calculations. This is because the observed temperature cannot trace the small temperature difference for a very short period in deep layers. This error can easily be corrected by estimating "true" temperature from running averages.
In order to estimate the length of the averaging window, an amplitude spectrum of the daily mean temperature at each level is calculated. If the spectrum at some frequency is smaller than the resolution, then the fluctuation at that frequency has no significance. Temperature resolution improves with the length of average window as:
where δT0 is the daily temperature accuracy, taken as 0.1 K, and N is the number of data to be averaged. An example of amplitude spectra of snow temperature is given in Figure 2, together with the accuracy limit, Equation 3. The spectrum of temperature intersects the accuracy limit at a frequency of about 0.08 day-1, which correspond to a period of about 12 days. The fluctuations in higher frequencies than this must be neglected because they cannot be measured in principle.
Fig. 2. Spectrum of temperature at the depth of 5 m (solid line) in comparison with the limit of accuracy when averaged over a period.
table 2. Length Of Running Average Window
In this case, the length of running average window is taken to be 15 days.
Table 2 lists the lengths of running average window at levels below 0.5 m which are determined by the method described above. The lengths are longer at deeper levels because only the lower frequency fluctuations must be present there.
Correction for the Integration Error
Another significant cause of the errors is the integration error b). In order to eliminate this error, the thermal difrusivity, K, for each layer is assumed to be constant and a better profile than the straight line is given by solving the well known diffusion equation ?T/?t = K(?2 T/?z2) (eg Munn 1966). Further, because of the linearity in the diffusion equation, any fluctuation in snow temperature can be expressed as a composition of sinusoidal fluctuations of various frequencies. If a sinusoidal temperature fluctuation, T = A0 cos ωt, is imposed at the interface z = 0, the temperature at depth z is given by:
where ω is the radian frequency, A_ is the amplitude of the sinusoidal temperature wave at z = 0 and 1 is defined by 
.
A contribution to the total integral from a layer of thickness, d, can now be defined for the straight-line profile as:
and for the constant diffusivity profile as:
The integral error b) can be corrected by comparing Is with Ic.
The contribution by straight-line profile, Is, can be calculated with
calculated with ?T/?t = -A0ωe(-Z/1)sin(ωt-z/1) and with a little trigonomic algebra as
where:
On the other hand, if the integral is calculated analytically, the contribution from the layer is calculated as
where
The integral formulae:
are used during the calculation. Therefore, an amplitude ratio RA≡|Ic|/IB| and a phase difference:
become, after a little calculation
Figure 3 displays RA and θ as a function of β=d/l by solid curves. The abcissa also shows ω/ω0, = ß2 where ω0 = 2K/d2.
Fig. 3. Amplitude and phase coefficients for correcting the integral error.
Calculation of RA and θ values by Equation 7 takes a long time on a microcomputer used in the present study. Therefore, Equation 7 is approximated as
which are shown by dotted lines in Figure 3. The asymptotic behavior of e by Equation 8 at small ω/ω0 deviates a little from Equation 7, because the coefficient (0.12) is chosen by fitting the curve in the intermediate region (1<ω/ω?<10).
Correction of the integral error is made as follows. First, a time series of Is is caiculated from the original data. Second, amplitude and phase spectra of I are calculated by using the fast Fourier transform (FFT) method. Third, the amplitude and phase at a particular frequency are corrected according to Equation 8. And at last, the time series of Ic is composed by inverse FFT from the corrected amplitude and phase spectra.
Calculation of Hourly Heat Flux Values
Hourly values of snow heat flux at the surface are calculated to see the validity of the corrections described above. Fast oscillations such as daily variations of temperature do not penetrate into layers deeper than, say, 0.5 m; thus the calculation is divided into two parts. First, daily heat flux at the depth of 0.5 m is calculated applying the correction on the temperature data with the running average method described above. Next, the hourly heat flux difference between the surface and 0.5 m depth is calculated with corrections described above. Adding the hourly values to the daily values at 0.5 m gives the heat flux at the surface.
The observation period is 312 days, thus giving 7488 hours of data. Because of lack of memory in the computer, the correction of the integral error is done by dividing the data into 8 blocks of 936 data each. In order to calculate FFT and inverse FFT, actual block length is taken at 1024, thus the blocks overlap with adjacent blocks.
Fig. 4. Hourly values of snow heat flux calculated with (lower) and without (upper) correction. The abscissa is day of the month (March 1979). The net radiation heat flux is shown by a dashed line for comparison. The heat fluxes are positive downward.
Figure 4 shows an example of calculated heat flux without corrections in the upper part and that with corrections by a solid line in the lower part. The results without corrections are filled entirely with noise and the effect of correction is remarkable.
Figure 4 also shows net radiation heat flux (Yamanouchi and others 1981) by a dashed line. Snow heat flux is very small compared with the net radiation but they agree in phase, which indicates that the corrections are reasonable. The remaining error in S0 is about 30%, from Figure I, but it is sufficient for the consideration of heat balance.
Concluding Remarks
There are other methods to estimate short term snow heat flux from the temperature measurement, including the one which uses the Fourier series of temperature (eg Dairymple and others 1966). In this case, the problem with integration, (Errors, b) above), does not occur but improper spacing of sensors will cause inaccurate estimates of the phase differences. Also, there still exists the error from factor a), thus the estimation of short term flux is difficult.
Another method would be to assume the thermal conductivity and to calculate S0 from the flux-gradient relations. But there is considerable scatter in the results on conductivity obtained by many workers (Mellor 1964). The present correction method also uses the thermal diffusivity K, which is the thermal conductivity divided by the density, but the correction is made on only a part of the data. Therefore, the errors in K do not affect the results seriously.
Consideration of the errors lead to a design principle of snow temperature observation. For example, the characteristic length 1 provides a measure of spacing of sensors. Or, a desired accuracy in the temperature derivative demands an accuracy in the temperature. But, an observation following the principle strictly may sometimes be expensive. In such cases, the methods presented here can be used to correct the errors.
Correction of errors c) and d) is not discussed here because their significance is smaller than the others. Repeated measurements of snow density will decrease the errors from c), whereas the method used by Dairymple and others (1966), who assumed genuine heat conduction in deep layer (z > 7 m), can be used to minimize the error from d).
The error and correction described above can also be applied to other field of conductive heat flux studies as long as the phase change of the snow is negligible such as in dry snow. If the snow depth is small and the integral can not be done to a very deep level where the heat flux vanishes, then the heat flux at the bottom or at some depth must be determined by another method such as using a heat-flux plate,
Acknowledgements
The authors thank the wintering party at Mizuho Station, led by Dr S Mae, and earlier discussions from Dr. N Ishikawa of the Institute of Low Temperature Science, Hokkaido University and Dr S Kawaguchi of National Institute of Polar Research.
