Modelling temperature variations in polar snow using DAISY

Abstract

A physics-based snow model has been calibrated using data collected at Halley Bay, Antarctica, during the International Geophysical Year. Variations in snow temperature and density are well-simulated using values for the model parameters within the range reported from other polar field experiments. The effect of uncertainty in the parameter values on the accuracy of the predictions is no greater than the effect of instrumental error in the input data. Thus, this model can be used with parameters determined a priori rather than by optimization. The model has been validated using an independent data set from Halley Bay and then used to estimate 10 m temperatures on the Antarctic Peninsula plateau over the last half-century.

Introduction

The surface snow cover in Antarctica forms a layer oforder 10 – 100 m deep over most of the continent. The density of the snow increases with depth to a value of about 800 kg m⁻³ at which point the air is contained in isolated pores. This pore cut-off density marks the transition to ice which can be up to 4000 m thick. It may at first sight appear that the relatively thin surface layer can be of no very great interest to Antarctic scientists. However, processes in polar snow have become important in several key areas of research:

• (i) Mass, energy and momentum transfers at the air snow interface are required as lower boundary conditions by atmospheric modellers working at all scales. Micrometeorological boundary-layer studies, investigations of katabatic winds, regional and global climate models all need to consider processes at least in the top few metres of the snow to derive accurate boundary conditions.

• (ii) Over much of Antarctica the only data collected on the ground have been density profiles of the top 10 m or so of the snow and the snow temperature at or near 10 m. These data have been used to construct climate maps of Antarctica but the derivation of surface air temperature from measurements of snow temperature depends on a correct understanding of heat flow in the snow cover.

• (iii) Ice-core data have been used very successfully to reconstruct atmospheric conditions in the past. The transfer functions relating to atmospheric and ice-core parameters are influenced by processes in the deposited snow as it is gradually transformed to ice.

• (iv) Remotely sensed passive microwave data can be used to determine surface temperature and accumulation provided that the density and grain-size for the first metres of snow are known.

Physics-based models for simulating processes in snow have been developed by many authors; those few which have been translated into fully developed computer programs include SNTHERM (Jordan, 1991), CROCUS (Brun and others, 1989, 1992) and DAISY (Bader and Weilenmann, 1992). In this paper, the performance of the DAISY program is tested using data from Antarctica and then the model is used to demonstrate the response of 10 m temperatures to climate change.

2. Data from Halley Bay Station

The Royal Society’s expedition to Halley Bay, Antarctica, was undertaken as part of the United Kingdom’s contribution to the International Geophysical Year (IGY). The measurements made at Halley Bay from 1956 to 1958 form an excellent data set for testing the capacity of snow models such as DAISY to simulate heat and mass transfer in polar conditions. They are all the more valuable because over the last 30 years a surface and upper-air observing programme has been maintained at Halley by the Falkland Islands Dependencies Survey (FIDS), later to become the British Antarctic Survey (BAS). This means that the meteorological conditions observed during the IGY can be analysed in the context of present-day understanding of the climatology of the area. Furthermore, comprehensive micrometeorological experiments in 1986 (King, 1994) and 1991 (King and Anderson, 1994) allow independent estimation of two key model parameters, aerodynamic roughness length for turbulent transfer and albedo, and validation of the model using an independent data set.

Synoptic Meteorological Data

Halley Bay station was established on the Brunt Ice Shelf at 75 °31′S, 26 °37′ W about 2 km from the ice front. Standard synoptic meteorological observations at 3 h intervals began on 8 January 1957 and were continued until 31 December 1958. The data are recorded in volume IV of the Royal Society International Geophysical Year Expedition Reports (table 237, p. 166) (MacDowall and others, 1964). The measurements used in the DAISY model are air temperature, humidity and wind speed.

From 1956 to 15 March 1957, air temperature and humidity were measured in a Stevenson screen located 45 m from the main hut on a bearing of 020 ° (Fig. 1, point A). The screen was then moved to a site 55 m from the main hut on a bearing of 080 ° (point B), and 9 m north of the anemometer tower (point C). From 1 February 1957 to 20 January 1958, the Assman Psychrometer used as the standard instrument for air temperature and humidity was mounted in the Stevenson screen at the standard height of 1.5 m above the surface. It was then moved to the base of the anemometer tower and maintained at 1.4 m above the snow surface. MacDowall and others (1964) pointed out that at air temperatures of lower than about −15 ° C the ice-bulb temperature (which can only be measured to a precision of about ±0.2 ° C) was not an accurate measure of relative humidity. However, we are concerned with the absolute magnitude of water-vapour and latent-heat transfers between air and snow. When temperatures are very low, the amount of water vapour in the air is negligible, whatever the relative humidity, so the Psychrometer errors do not have a significant effect on the simulation of the mass and energy balance at the snow surface. The model results discussed later (in section 4) show, for example, that in the winter fluctuations in air temperature from −20 ° to −50 ° C are associated with less than 1 W m⁻² change in the latent-heat flux; in the summer, air-temperature fluctuations from 0 ° to −10 ° C are associated with latent-heat flux changes of the order of 20 W m⁻²

Fig. 1. Sites of measurements made at Halley Bay during the IGY. A and B, sites of Stevenson screen; C, anemometer mast; 1 and 2, accumulation stakes read monthly; a – f, accumulation stakes read daily.

Continuous records of wind speed began on 6 March 1957 when a cup anemometer linked to a recording voltmeter was installed on the tower (point C in Figure. 1). The values reported at the synoptic hours were 10 min averages for a nominal height of 10 m. The observers removed hoarfrost deposits from the cups when necessary, so these observations are probably reasonably accurate, except perhaps at low wind speeds after a storm when snow in the gears might hinder rotation.

Radiation

Incoming solar radiation was measured using a solarimeter mounted on the roof of the main hut (Fig. 1) from August 1956. Net all-wave radiation was measured from March 1957 using a ventilated flux-plate radiometer mounted about 1 m above the snow surface between the hut and anemometer tower (C). The data are recorded in volume III of the Royal Society International Geophysical Year Expedition Reports (MacDowall and Tribble, 1964, p. 123 – 59). MacDowall and Tribble discussed the expected errors in these radiation measurements in some detail. The daily incoming solar radiation (used in DAISY) is estimated to be accurate to 4% + (2 × 10⁵)J m⁻² and the daily net radiation to 15% + (4 × 10⁵)J m⁻², although this latter figure may be a little pessimistic.

Snow Density

Snow density was measured by taking samples of known volume from the side walls of pits. Figure. 2 shows profiles measured in May 1957 and December 1958. The May 1957 profile begins at the 1955 – 56 summer surface (i.e. around January 1956); densities from the top layer covering the period January 1956 – May 1957 were not recorded. Summer layers produced by melting can be distinguished and have been used by MacDowall (1964) to identify the annual layers. He considered, for example, that 120 cm of snow lies between the January 1957 and January 1958 summer layers. The error in the density measurements is of the order of ± 50 kg m⁻³.

Fig. 2. Snow-density profiles (after MacDowall, 1964).

Snow Temperatures

Temperatures in the snow were measured using thermocouples. On 1 May 1957 seven thermocouples were placed at depths of 0.16, 0.34, 0.60, 1.21, 3.05, 6.10 and 12.19 m. MacDowall (personal communication, 1994) has confirmed that cables ran directly from the Main Hut (Fig. 1) to the thermocouples which were free to move downwards at the same rate as the surrounding snow.

MacDowall (1964) interpolated and smoothed the temperature records produced by the seven thermocouples to produce temperature curves for three fixed depths below the surface of the snow. In order to do this, he tacitly assumed that the thermocouples all moved downwards as the snow compacted in exactly the same way as the stake used for the level measurements. However, if there was any differential movement between the thermocouples and the stake, the curves shown in his Figure. 13 are not valid. The most likely scenario is that each thermocouple moved downwards at the compaction rate of the layer of snow where it was inserted, whereas the measuring stake moved downwards at the (slower) compaction rate at its base.

Accumulation

During the period 15 June 1956 – 31 December 1958, accumulation measurements were made at three separate sites (Fig. 1); measurements were made at irregular intervals at site 1, monthly at site 2 and at daily intervals at sites a – f. First runs with DAISY were made using the data from site 2, where the accumulation over the year from May 1957 was 100.5 cm of snow. This is comparable with the 120 cm between January 1957 and January 1958 estimated from the snow pit.

New Snow Density

The density of new (i.e. dry, falling) snow was not measured directly in the IGY experiment but an estimate of nearly-new snow density (after initial settling) may be made by examining the densities measured at the accumulation site when the mean height change in a day is recorded as more than 1 cm. The mean of 32 values over the period is (320 ± 70) kg m⁻³ with the minimum value being 200 kg m⁻³.

3. Daisy

DAISY has been described in detail in Bader and Weilenmann (1992), where the complete set of equations defining the model was given. One of these is the heat-flow equation, which governs variations in temperature, T, and freezing rate, m¹², with time, t, and height above the ground, z:

(1)

Here, superscripts (1) and (2) refer to ice and water respectively, ρ^(w) is the partial density, v^(w) is the velocity and C_p ^(w) is the specific heat of component w. The total density is

where ρ⁽³⁾ is the density of dry air and the effective heating Capacity is similarly defined as

L₁₂ is the latent heat of freezing. The water content is assumed to be zero if T < 0 ° C.

The heat-flow equation includes an effective thermal conductivity, λ_eff which depends on the density of the snow. Mellor (1977) gave the variation of experimentally determined values of effective conductivity with the density shown in Figure. 3 (shaded area). Morris (1983) has proposed the expression

(2)

shown in Figure. 3 as the dashed line. Although this expression gives the broad pattern of the increase in effective thermal conductivity with density, it is clear that the uncertainty is large. The DAISY model allows λ_eff to be adjusted by multiplying by a factor K₁ Curves for K₁ = 0.5 and 2.0 are shown by the solid lines in Figure. 3. MacDowall (1964) calculated thermal diffusivity from the temperature curves shown in his Figure. 13 and found values of 6 × 10⁻⁷ m² s⁻¹ at the surface to 12 × 10⁻⁷ m² s⁻¹ from 3 to 6 m down. Equation (2) gives values of 5.5 × 10⁻⁷ m² s⁻¹ for ρ = 500 kg m⁻³ and 4.4 × 10⁻⁷ m² s⁻¹ for ρ = 300 kg m⁻³ , so MacDowall’s field data suggest K₁ should be between 1 and 2. However, we have already suggested (in section 2) that MacDowall’s temperature curves are based on a false assumption, so his calculations of diffusivity should be treated with caution. Anderson (1994) has calculated a thermal diffusivity of 5 × 10⁻⁷ m² s⁻¹ for the upper 1.5 m of the snow cover at Halley between March and September 1991.

Fig. 3. The effective thermal conductivity of snow (after Mellor, 1977).

Equation (1) also includes a source term, R_s, which arises from the absorption of solar radiation within the snow cover. R_s is calculated using an extinction coefficient β which is expected to vary within the range 20 – 160 m⁻¹. β increases with density and decreases with increasing grain-size (Mellor, 1977).

The source term R_w in Equation (1) is included to account for energy produced by wind-pumping (Colbeck, 1989) and depends on the wavelength and height of sastrugi. It is set to zеrо for this modelling exercise.

The program calculates the change in density of each layer with time using the densification equation suggested by Bader (1960, 1962)

(3)

where ϵ is the strain rate, σ is the overburden pressure and η is the compactive viscosity. This is, of course, a highly simplified representation of the complex viscoelastic behaviour of snow but Equation (3) is an appropriate first formulation given the other simplifications in the model. It is assumed that the compactive viscosity of a layer varies with its density and absolute temperature according to the expression

(4)

where H₀, C and E (the activation energy) are constants. R is the gas constant (8.314 J mol⁻¹ Κ⁻¹). Figure. 4 shows the range of values of η obtained in field experiments on natural snow packs (Mellor, 1975). In general, the data are consistent with Equation (4). Mellor remarked that the discrepancy between values for polar snowfields and seasonal snowfields seems to be too great to be explained by temperature difference alone and suggested that the seasonal snow values reflect basic shortcomings in the concept that densification occurs by slow viscous creep. In his view, densification occurs more by a quasiplastic collapse over a limited time after each significant increase in stress. Thus, apparent values of η are partly dependent on the interval between snowfalls. If large discrete snowfalls occur at frequent intervals on a seasonal snow-pack, the apparent value of η will be lower than expected. However, there is also a difficulty in defining the temperature at which the densification process is taking place in these field experiments. The new, low-density snow at the surface of the snow-pack is subject to a fluctuating temperature on the diurnal time-scale. The effective temperature for the densification process, i.e. the temperature to be used in Equation (4), will be higher than the average temperature over periods longer than a day and certainly higher than the mean annual average temperature which at most polar sites is all that is known.

Fig. 4. The compactive viscosity of snow as a function of density (after Mellor, 1975). The shaded areas show experimental data from; (A) Greenland and Antarctica, −20 ° to −50 ° C (Bader, 1953); (B) Seasonal snow, Japan, 0 ° to −10 ° C (Kojima, 1967); (C) Alps and Rocky Mountains (Keeler, 1969); (D) Uniaxial strain-creep tests, −6 ° to −8 ° C (Keeler, 1969); (E) Uniaxial strain-creep tests, −23 ° and −48 ° C (Mellor and Hendrickson, 1965). The curves show the predicted variation using Equation (4) with parameter values given by Bader and Weilenmann (1992) for (a) T = 0 ° C, (b) T = −20 ° C, (c) T = −50 ° C, and by Kojima (1964) for (d) T = 0 ° C, (e) T = −20 ° C, (f) T = −50 ° C.

Bader and Weilenmann (1992) set the parameters in Equation (4) to the values H₀ = 0.18 × 10⁻⁵kg m⁻¹ s⁻¹, C = 0.02 m³ kg⁻¹ and E/R = 8110 K. which gives the variation of η with ρ at temperatures T = 0 °, −20 ° and −50 ° C shown in Figure. 4 (solid lines). The fit is good for the seasonal snow data set B but the predicted values of η are somewhat too low for the polar data set A. Kojima (1964) has studied the densification of snow in Antarctica in considerable detail and suggests an activation energy E = 12 k cal mol⁻¹ = 50.2 kJ mol⁻¹ for C = 0.024 m³ kg⁻¹ and H₀ = 5.38 × 10⁻³kg m⁻¹ s⁻¹, hence E/R = 6042 K. These values give curves (dashed lines) which are a better fit for the polar data set A but the predicted values of η are too high for the seasonal snow-pack. In this application of DAISY, the Kojima parameter values are used in Equation (4) and η is further adjusted by multiplying by a factor K₂ .

Mass transfer by phase change and the transport of water vapour within the snow-pack are not dealt with explicitly within the model. The use of an effective thermal conductivity and empirical densification equation will allow these processes to be accounted for to a certain extent but it must be accepted that DAISY, in its current form, cannot simulate some phenomena linked to vapour transport, such as the formation of internal layers of hoar frost.

Boundary Conditions

The upper boundary condition consists of an equation for the energy flux, ε_z, across the air-snow surface. This is calculated as the sum of net longwave radiation, sensible-heat transfer, latent-heat transfer and heat convected by precipitation. Where net all-wave radiation, R_n, and incoming solar radiation, S↓ are available, as at Halley Bay, the net longwave radiation is calculated as

where α is the albedo of the snow surface and may be expected to vary from 0.75 for old polar snow in summer to 0.95 for fresh dry snow (Mellor, 1977). Gardiner and Shanklin (1989) gave mean hourly values over the period 1963 – 82 at Halley ranging from 0.79 to 0.84.

The equations for turbulent transfer of sensible and latent heat used in DAISY are functions of the aerodynamic roughness length z₀ and the roughness lengths for heat-and water-vapour flux z_T, z_q. King (1994) has calculated z₀ from measurements made at 5 m above the snow surface at Halley in near-neutral conditions. His value of z₀ = (1.1 ±0.1) × 10⁻⁴ m is consistent with other measurements over flat snow surfaces (Morris, 1989). In a later experiment (STABLE II) at Halley, King and Anderson (1994) found a value of z₀ = 5.6 × 10⁻⁵ m with a range of uncertainty from 5.0 × 10⁻⁵ m to 6.1 × 10⁻⁵ m and a value of z_T = 5.0 × 10⁻⁵ m with bounds of 1.9 × 10⁻² m and 1.3 × 10⁻¹ m. There were insufficient data to determine z_q but the measurements were not incompatible with the generally accepted hypothesis that z_T ≈ z_q.

The turbulent-transfer coefficients are adjusted for the stability of the boundary layer of the atmosphere by multiplying by a factor

where Ri is the Richardson number calculated at the standard height of 1.5 m, at which air temperature and humidity are measured. The stability correction factor is limited to the range 0.84 ≤ K₃ ≤ 100 following Brun and others (1989), who proposed an empirical linear equation for K₃ as a function of wind speed. The minimum value of K₃ varies from site to site but a typical value is 0.8 (Bader and Weilenmann, 1992). There is clearly scope for investigating alternative expressions for the stability corrections and this will be pursued in a further paper. Since the wind speed is measured at 10 m the value at 1.5 m needed to calculate Ri is estimated by assuming that the wind-speed profile is logarithmic. In non-neutral conditions, this introduces some error. The model could certainly be improved by iteration to correct the estimate of wind speed given an estimate of Ri.

It is assumed that the temperature of the precipitation (in this case always snow) is the same as the air temperature. The initial density of a new snow layer is set to a fixed value which is expected to lie in the range 300 – 400 kg m⁻³.

4. Simulation of IGY Data

Initial Conditions

The calibration period was chosen to be 1 year starting from 1 May 1957. The initial density profile was calculated on the basis of the data from the two snow pits dug in May 1957 and December 1958. In December 1958, the thickness of the layer of snow which fell between January 1956 and May 1957 was measured as 120 cm. We assume therefore that in May 1957 the unrecorded upper part of the pit was 120 cm thick. The density profile in this surface layer can be reconstructed using the data from the top 120 cm of the December 1958 pit as a guide. In the top 50 cm, we assume that the May values will be 30 kg m⁻³ less than the December values to take account of compaction. Table. 1 shows the measured and estimated densities at 20 cm intervals.

Table 1. Initial structure of the snow cover

The next step is to define the layers which may be considered to have the same density. The maximum number of layers (including those added by accumulation during the year) is around 35 for a PC with 4 M byte storage so that the initial number of layers was chosen to be 18. More detail is required in the upper part of the Snow-Pack where temperature variations are larger. Table. 1 shows the layers chosen and the average density for each layer. The spread of densities in one layer is not allowed to be more than 80 kg m⁻³ (the measurement error in one density measurement is ±50 kg m⁻³). The number of grid points in each layer must be at least 2 and is chosen to produce a spacing which varies from 2 to 40 cm, according to the density and temperature gradients expected in the layer.

It should be noted that, since some data from the December 1958 snow pit have been used as a guide to reconstruction of the initial conditions, the measured density profile cannot be regarded as a completely independent criterion for validating the simulation of density when DAISY is run on to December 1958. However, the measured density profile in the top layer of snow which fell after 1 May 1957 can of course be used to check the simulated densities, since these will be independent of initial conditions.

The initial temperature profile was derived from the profile measured at 1800 h on 1 May 1957 after the newly installed thermocouples had had time to adjust to the temperature of the snow. The measured temperatures were adjusted so that they conformed with the constant-temperature lower-boundary condition, T = −18.4 ° C at 12 m depth.

Input data

Early runs with DAISY were made using the accumulation rates measured at the stake farm. It soon became clear that these data were not applicable to the thermocouple site and indeed MacDowall (1964) remarked that at the thermocouple site “the rate of accumulation was greater than that normal over the level ice-shelf owing to local drift effects”. It therefore was necessary to try to reconstruct the accumulation record at the thermocouple site using the measurements of the level of the snow surface made at that site plus an estimate of the average surface-settling rate. The snow surface was measured relative to a reference mark on a Dexion mast, which may itself have been sinking into the snow. The choice of input-accumulation curve for a given set of parameters used in DAISY is made so that the simulated snow-surface level agrees as closely as possible with the level measured with respect to the Dexion mast.

Results

Figure. 5 shows results from DAISY run HY with albedo α = 0.9, extinction depth β = 20 m⁻¹, aerodynamic roughness length z₀ = 10⁻⁴ m, roughness-length ratio z_T/z₀ = z_q/z₀ = 1, new-snow density 400 kg m⁻³, compactive viscosity given by Kojima’s parameters in Equation (4) for η and effective thermal conductivity given by Equation (2) for λ_eff multiplied by 1.8. The initial snow-density profile for this calibration run was taken to be that given in Table. 1, except that, since the density of a new layer is taken to be 400 kg m⁻³, the initial densities of layers 17 and 18 were increased to 400 kg m⁻³ to be compatible with this assumption.

Fig. 5. Simulations produced by DAISY run HY with parameters as shown in Table. 3. (a) Temperatures at fixed depths below the snow surface; (b) Temperatures at fixed layers. The measured temperatures from thermocouples in the layers are also shown; (c) Height of the snow surface and internal layers above the datum; (d) Density changes in the top metre of the snow-pack.

Figure. 5a shows the simulated temperatures at fixed depths of −1.5, 3.0. 4.5 and 6.0 m below the surface. Comparison with the temperature curves shown in MacDowall’s Figure. 13 (MacDowall, 1964) shows that the two sets of curves are similar for the first part of the year. After November, MacDowall’s temperatures are lower than the DAISY temperatures. We believe this is because MacDowall’s calculations were based on the erroneous assumption that the thermocouples were not moving downwards with respect to the measuring stake. Figure. 5b shows the temperatures for fixed layers compared to the measured temperatures for the thermocouples installed in each layer. The fit is good except for a period in August and September when the top three thermocouples (2,3 and 4) recorded larger temperature variations than predicted. We believe that, for most of the year, the thermocouples moved with the snow layers. The fact that 2,3 and 4 seem to be closer to the surface than predicted in early spring may indicate that these thermocouples were restricted by their cables for a while or that the compaction rate over the winter-accumulation period has been overestimated.

Figure. 5c shows the height of the snow surface measured at the Dexion mast and simulated by DAlSY. The position of the measured points depends on the assumption that the Dexion mast was fixed with respect to the layer at 12 m depth. The initial positions of the thermocouples are also shown. Since the measured heights have been used to derive accumulation as an input to DAISY, the two surface curves are not independent where their gradient is positive. In summer, however, during periods without snowfall, the curves are independent and the close fit between measured and simulated surfaces suggests that surface settling and compaction are correctly simulated. The figure also shows the evolution of layers in the snow cover. The assumption that the layer at 12 m depth can be regarded as a fixed height datum is clearly reasonable but it is also clear that, if the Dexion mast had only been inserted 2 m into the snow at the beginning of the period, it would have moved about 1 m downwards after 1 year. Unfortunately, the dimensions and method of erection of the mast have not been recorded. By the end of the year, the increase in surface height above the datum at 12 m depth is 140 cm. (Thus, if the ice shelf is in equilibrium, the ice velocity at the datum level should have a downward component of about 1 m year⁻¹, so the height of the ice-shelf surface relative to sea level is maintained.) This is compatible with the local accumulation rate of about 1 m year⁻¹. However, by May 1958, there is in fact a 3 m layer of snow above the initial snow surface of May 1957. It could be that this high accumulation rate is genuine (the site was susceptible to drifting) but it is also possible that the Dexion mast was moving downwards. This would mean that the accumulation has been overestimated.

The density profile in the top 3 m layer is independent of the initial conditions input to DAISY and can be used as another check on the compaction rate. The density changes at fixed depths below the snow surface are shown in Figure. 5d. After a new snowfall, a point at a given depth below the surface lies in a new layer, so the density falls abruptly. Then, as the layer compacts, the density gradually rises. In the period from July to September there is a large amount of snowfall, so the density of the upper 1 m is mostly controlled by the choice of the new-snow density. Later in the year, the choice of compaction rate has more influence. By the end of 1957, the new snow has developed a density profile very similar to that observed in the top 1 m of the December 1958 snow pit, good evidence that the compaction rate used in run HY is correct. However, DAISY does not reproduce the ice layers which are shown in Figure. 2 because in the current version of the model meltwater is not allowed to move downwards through the snow and then refreeze. Ice layers are not included in the initial density profile and do not form in the new snow.

Table. 2 shows the components of the energy balance at the upper boundary averaged over each month calculated by DAISY run HY. Energy flux from the air to the snow is positive. The figures in brackets are seasonal means measured in 1991 (King and others, in press), or in the case of net radiation, the 1963 – 82 mean (Gardiner and Shanklin, 1989), measured at the same site and quoted for comparison.

Table 2. Mean monthly components of the energy balance at the upper boundary

In general, the values for net radiation and sensible-heat flux are comparable. However, in 1957, the monthly mean latent heat flux was always negative (heat was lost from the snow because of sublimation and evaporation), whereas in winter 1991 there was a small latent heat gain from condensation. Note that for 3 months in the summer DAISY calculates the latent-heat flux to be the largest component of the energy balance at the upper boundary. The mean annual energy flux at the upper boundary over the year is −10.7 W m⁻², which counteracts the mean rate of absorption of solar radiation in the snow, which is 10.5 W m⁻².

5. Sensitivity Analysis

The main advantage of a physics-based model is that it is possible to set limits on the expected range of the parameters and distinguish the individual effects of each one. Figure. 6 shows the envelope of temperatures at a fixed depth of 1.5 m below the surface produced by varying one parameter at a time about the values used in the calibration run, that is, by the runs shown in Table. 3.

Fig. 6. Sensitivity of temperature at 1.5 m below the surface to changes in effective thermal conductivity, extinction coefficient, albedo, aerodynamic roughness length, new snow density and compactive viscosity.

Table 3. Parameter values used in the sensitivity analysis

Figure. 6 shows that the temperatures at 1.5 m given by runs HP, HZ, HF, HB and HC are never more than 1 ° C from the temperature given by the calibration run HY. The uncertainty is about 10% of the amplitude of the annual wave. This value is, in fact, somewhat less than the uncertainty expected, because of errors in input data (especially net radiation, which is only measured to 15% accuracy). Thus DAISY’s performance in simulating temperature is not limited by uncertainty in parameter values.

Decreasing the effective thermal conductivity decreases the amplitude of the temperature wave at all depths and increases the lag time. For this reason, run HZ shows most deviation from HY, since the temperature curve is decreased in both amplitude and offset. Since net radiation is measured, the choice of albedo and extinction coefficient does not affect the magnitude of the radiation input, merely where it is absorbed. The reduction of albedo to 0.8 in run HA means that more of the net radiation is taken to be shortwave and hence the temperature at 1.5 m is greater in the austral summer by about 1 °С. Increasing the extinction depth for shortwave radiation (run HB) also increases the temperature in the austral summer by about 0.5 ° C.

Because the sensible-heat flux is quite small, the temperature of the snow surface is close to the air temperature for most of the year and thus the simulations are not very sensitive to variations in roughness-length ratio. Increasing the ratio to 500 (run HC) produces a maximum reduction in temperature at 1.5 m of about 1 ° C in January when the snow-surface temperature is warm because of the positive net radiation (see Table. 2) but the air temperature is about 2.5 ° C lower and the boundary layer becomes unstable. Decreasing the new-snow density (run HP) leads to warmer temperatures at 1.5 m in July and August, again by about 1 ° C. This is because the thermal conductivity of the overlying snow is lower and the winter cold wave is attenuated. In the summer, there is less snowfall so the effect of the choice of new-snow density is less significant. Decreasing the compactive viscosity by a factor of 10 (run HF) has negligible effect on the temperature at 1.5 m. This is because, in order to maintain the correct snow surface height the surface settling rate used in run HY (7.8 × 10⁻⁸ m s⁻¹) is reduced to 1.3 × 10⁻⁸ m s⁻¹ in run HF. Thus, the density changes in the top 1.5 m of snow are relatively small and the thermal conductivity remains much the same. Because reliable independent measurements of accumulation are not available from the thermocouple site, it is difficult to optimize the compactive viscosity parameter. However, it is possible to say that K₂ cannot be less than about 0.1, since a negative surface-settling rate would then be required to maintain the correct snow-surface height.

Comparison of final density profiles for runs HY, HP and HF shows that reducing either the new-snow density or the compactive viscosity reduces the predicted density in the upper 3 m of snow, though not by more than 50 kg m⁻³, which is the minimum expected error in the measured initial densities. Run HF gives a good fit to the profile below 2 m, whereas the calibration run HY gives densities which are slightly too high, though again not by more than 50 kg m⁻³. Thus, DAISY’s performance in simulating density is not limited by uncertainty in parameter values.

Table. 4 shows the annual average components of the energy balance at the upper and lower boundaries and for the snow-pack as a whole. Runs HZ and HA were terminated because of numerical instabilities after 355 d and run HY after 364 d. All other runs were for the full 365 d. For this reason, the calculated values of annual average net longwave radiation for runs HY and HZ are slightly different (0.1 and 0.2 W m⁻² less negative) than the value of −11.3 W m⁻² for the other runs with the same albedo. For all runs, except HC, the turbulent transfer of sensible heat into the snow is nearly balanced by turbulent transfer of latent heat of sublimation to the air. The shortwave radiation absorbed by the snow is balanced by longwave radiation to the air. Heat added to the snow by precipitation and the heat flux across the lower boundary are small. Run HC, with roughness-length ratio z_T/z₀ = 500, produced a large decrease in the heat content of the snow cover over the annual cycle because of the increased evaporation. This is very much an extreme case and we suspect that further work at Halley, now being undertaken (personal communication from J.C. King), may lead to downward revision of his estimate for the roughness-length ratio at this site.

Table 4. Parameter values used in the sensitivity analysis

6. Validation

This paper has described calibration of the DAISY model using IGY data from Halley Bay. In order to validate the model, we have used a completely independent data set from Halley collected during the STABLE II experiment in 1991. This work is described in detail by Morris and others (1994); here, it is only necessary to remark that by using the “effective” values of the DAISY parameters, i.e. those used in the calibration run HY, good simulations of the STABLE II data were obtained, except during short periods of very high stability in the atmospheric boundary layer. Thus we are encouraged to believe that DAISY can be used for prediction with parameters determined a priori rather than by optimization.

7. Daisy in Use: an Example

To demonstrate the potential of physically based snow models, such as DAISY, we describe here briefly how the model has been used to estimate the variation of 10 m temperature on the Antarctic Peninsula plateau over the 49 years for which meteorological data from Faraday station are available. During this period, there has been a systematic longterm warming of 2 ° C at Faraday (Stark, 1994) but also high inter-annual variability in air temperature (King, 1994). Measurements of 10 m temperature in boreholes have been made by glaciologists since 1960 and used as an indication of the mean annual surface temperature. By using DAISY we can assess the significance of the sparse historical data from the Antarctic Peninsula plateau and whether, in general, 10 m temperature measurements can be used to detect climate change in this region where meteorological data have not been collected in the past.

The upper-boundary condition at a given site on the plateau is chosen to be

where [] is the mean monthly air temperature at Faraday measured at a fixed height, nominally 1.5 m, above the snow surface and Т₀(ϕ, λ, Η) is an offset temperature determined for each site from a lapse-rate equation for western Antarctic Peninsula sites north of 74 ° S (Morris and Vaughan, 1994). ϕ, λ and H are the latitude, longitude and altitude of the site. We choose a Dirichlet upper-boundary condition (specifying T_s ) rather than a Neumann condition (specifying (∂T/∂z)_s,) because the wind speed, humidity and net radiation needed to calculate the temperature gradient at the snow surface on the plateau cannot reasonably be estimated from the Faraday measurements. The source term R_s in Equation (1) is set to zero. All net solar radiation is assumed to be absorbed in the surface layer of snow rather than the top metre or so. This assumption will have little effect on the temperature variation at 10 m depth.

For the plateau sites, the mean annual temperature is around −20 ° C so the surface temperature T_s(ϕ, λ, H, t) will only rarely reach 0 ° C. It is not necessary to simulate water movement in the snow, although melting and refreezing are included in the model. Since temperature changes over 49 years will penetrate around seven times deeper than the annual wave, the lower boundary (at which the temperature is considered to be constant) is taken at a depth of 100 m. The initial density profile chosen for the simulation is given by Table. 5 and is based on that of an ice core from the region.

Table 5. Initial density profile

Accumulation measurements are very sparse over the 49 year period, so a constant average rate of 4 × 10⁻⁸ m (snow) s⁻¹, i.e. 0.517 m water a⁻¹, has been used in the simulation. The new-snow density was taken to be 410 kg m⁻¹. Figure. 7 shows the sum of the 10 m temperature, T₁₀, and the appropriate offset temperature as a function of time. The curve simulated by DAISY shows the annual cycle (attenuated from about ±11 ° C to ±0.3 ° C) superimposed on decadal variation of about ±1 °С lagged by 1 – 2 years with respect to the surface air temperature. So, for example, the minima in 10 m temperature in 1960 and 1988 appear in the surface record of annual mean temperature in 1959 and 1987 (King, 1994). The points show measured values from four areas on the plateau. 10 m temperature was measured at Charity Depot (site JB1) (70.0167 ° S, 64.4833 ° W, 2131 m), in 1974 by Bishop (Reynolds, 1981), in 1979 by Paren (Reynolds, 1981) and nearby, at Charity level line (69.9813 ° S, 64.2308 ° W, 1990 m), in 1992 by Morris (Morris and Mulvaney, 1966). These points are all labelled “C” in Figure. 7. In 1992, Morris also measured 10 m temperature at St Pancras (71.7430 ° S, 63.9963 ° W, 1861 m) (Morris and Mulvaney, 1996) close to the site JB3 (71.7000 ° S, 64.0833 ° W, 1886 m) where Bishop measured 10 m temperature in 1975 (Reynolds, 1981). These two points are labelled “P”. A less well-matched set of measurements are Morris’s 1992 temperature from Temnikow Nunataks (70.5673 ° S, 64.2443 ° W, 1600 m) (Morris and Mulvaney, 1996), Bishop’s 1975 temperature from site JB2 (70.8333 ° S, 64.4500 ° W, 1987 m) (Reynolds, 1981) and Mumford’s 1976 temperature from Bullihole (70.8833 ° S, 64,9500 ° W, 1835 m) (Reynolds, 1981), all labelled “T”. The earliest measurements available from the region come from the Antarctic Peninsula Traverse in 1962 (Shimizu, 1964). The closest modern site to Shimizu’s three traverse sites APT668 (73.900 ° S, 69.4333 ° W, 1215 m), APT700 (73.5333 ° S, 68.6167 ° W, 1045 m) and APT732 (73.7167 ° S, 67.2667 ° W, 1575 m) is Bishop’s site JB5 (73.700 ° S, 64.7833 ° W, 2007 m) (Reynolds, 1981). These data are all labelled “J” in Figure. 7.

Fig. 7. Snow temperature at 10 m depth on the Antarctic Peninsula plateau predicted by DAISY for the period 1 March 1944 – 1 March 1993 with measurements from sites near Charity Depot (C), Bishop site JB5 (J), Temnikow (T) and St Pancras (Ρ)

The two well-matched groups “C” and “P” follow the curve reasonably well as do two of the “T” group. However, the third “T” point, the 1975 value of (T₁₀ + T₀) from site JB2, is about 0.9 ° C lower than expected. The 1975 site is 1040 m higher than the 1992 site. Given the uncertainty in the altitudinal lapse rate of ±0.0007 ° C m⁻¹ (Morris and Vaughan, 1994), the estimate of the offset T₀ could be too low by 0.7 ° C. This is sufficient to explain why this site does not fit the curve. The same argument probably applies to the “J” group of points, where there is a marked difference in altitude between the APT sites and JB5. The fact that the three APT points are scattered is an indication of the uncertainty in all three spatial lapse rates; the measurements were made at the same time and, given correct values of T₀, should lie at the same point on Figure. 7. Furthermore, we must accept that the variation in air temperature at Faraday will be an imperfect representation of the variation in surface temperature on the plateau. The “J” points are farthest from Faraday, so least likely to follow the simulated 10 m curve.

Despite these caveats, the results are interesting as they demonstrate the importance of recognising that 10 m temperature measured at a given site is a function of the climate history and is not necessarily equal to the mean annual surface temperature for the same year. The variation simulated using the Faraday temperature record is large enough to explain the differences in measurements made at the same sites at different times and to give confidence that repeat measurements of 10 m temperature can be used to detect climate variation on the Antarctic Peninsula plateau. Further work, using temperature records derived from deep ice cores from the plateau and the complete set of 10 m temperature data available for the region, will be reported elsewhere.

8. Conclusions

The DAISY model has been calibrated for polar conditions using data from Halley Bay primarily because the quality of the meteorological data collected during the IGY is so good. Snow-pit measurements are available to initialise the model and test its capacity to simulate densification, and temperature was recorded at several levels in the snow throughout the IGY. The only weakness in the data set arises because the positions of the snow surface and thermocouples were not recorded with respect to a fixed level. Simulated snow temperatures at the estimated thermocouple positions agree well with the measured temperatures except for a 3 month period when it appears that the near-surface thermocouples are higher in the snow-pack than expected. Furthermore, the accumulation at the thermocouple site is much higher than the average local value, possibly as a result of drifting, but possibly because the snow surface was unwittingly measured with respect to a moving datum. However, these are relatively minor problems given the advantages of this extensive and well-documented collection of data.

The densification of the snow was well represented by an equation based on Kojima’s (1964) studies of Antarctic snow. Heat flow was best modelled assuming effective thermal conductivities in the upper part of the range reported by Mellor (1977). This probably reflects the relative importance of vapour diffusion and convective air and vapour movement in polar conditions. Good estimates of the energy input at the upper boundary of the snow could be made using an albedo of 0.9 and aerodynamic roughness length of 10⁻⁴ m. Thus all the parameters in the DAISY model have “effective” values which are as expected from independent studies. We have reported elsewhere (Morris and others, 1994) on the successful validation of the model using these effective values to simulate temperature variations during the 1991 STABLE II experiment at Halley.