Introduction

It is commonly assumed that the 10 m temperature of dry polar firn and the mean annual air temperature are simply related (Reference PatersonPatcrson, 1981, p. 188; Hookc and others, 1984). Observations at Agassiz Ice Cap, Ellesmere Island, Arctic Canada (Reference Clarke, Fisher and WaddingtonClarke and others, 1987; Reference Waddington, Magnusson, Fisher and KoernerWaddington and others, 1988,1989) reveal the existence of microclimate zones that cast doubt on this assumed simple connection. Near the crest of Agassiz Ice Cap, 10 m temperature varies systematically by up to 5°C over a distance of several kilometres. The variations appear to be related to wind patterns associated with the summit dome of the ice cap. The anomalous microclimate signal exceeds 50% of the generally accepted glacial-interglacial temperature difference in polar regions.

Wind pumping (Reference Clarke, Fisher and WaddingtonClarke and others, 1987) is one of several working hypotheses to explain the observed temperature variations at Agassiz Ice Cap. Air, a viscous fluid, can be forced through the firn, a rigid porous material, in response to surface air-pressure fluctuations associated with turbulence. The resulting frictional dissipation can heat the firn. Using a one-dimensional model, Reference Clarke, Fisher and WaddingtonClarke and others (1987) estimated that wind pumping could account for as much as 1–2°C of the thermal anomaly at Agassiz Ice Cap, if the magnitude and spatial variability of wind speed approach meteorologically extreme values. To assess the importance of frictional heating when air is forced to move through the rigid firn matrix, Reference Clarke, Fisher and WaddingtonClarke and others (1987) derived a theory for the penetration of one-dimensional (plane wave) pressure fluctuations into a permeable firn half-space (Fig. 1a). Subsequently, Col-beck (1989) derived equations for steady air flow in a snow layer having sinusoidal surface topography. In Col-beck’s analysis, topographic features on the snow surface generate surface air-pressure perturbations that are spatially and temporally fixed. This is a restricted example of two-dimensional wind pumping (Fig. 1b) in which both temporal and spatial variations in pressure would be permitted.

Fig. 1. Schematic diagram showing the distinction between one-, two- and three-dimensional wind pumping.

Recent measurements of pressure fluctuations above and within polar firn now suggest that the one-dimensional pressure wave assumed by Reference Clarke, Fisher and WaddingtonClarke and others (1987) is unrealistic at the frequencies (0.1–10 Hz) thought most likely to generate measurable frictional heating in the firn. Figure 2 shows measured power spectra and best-fitting straight-line approximations obtaincd from simultaneous pressure recordings at the snow surface and at 0.10 m depth in Agassiz Ice Cap. At the time of this determination the wind speed was 5.1 ms⁻¹at 0.60 m above the surface. In Figure 2 the broken line shows the power predicted at 0.10 m depth using the one-dimcnsional theory and measured values of permeability and density at the experimental site. Evidently, the actual attenuation with depth far exceeds that predicted by the one-dimensional theory. We attribute this discrepancy in part to horizontal inhomogeneity in the surface pressure field (Fig. 1c). Reference BergenBergen (1980) extended the standard von Karman logarithmic velocity-profile theory (e.g. Reference TrittonTritton, 1988, p. 342) to the case of a permeable boundary. Bergen’s table 1 provides data that allow us to estimate mean wind speeds of 0.4–0.8 m s⁻¹ at the surface of a snowpack when his measured wind profile showed a few metres per second at 0.60 m height, comparable to our measured wind speed at Agassiz Ice Cap. Spatial variability of pressure cannot be measured without an array of sensors. However, for advection-dominated flows, the stochastic properties of the spatial variability at a given instant can often be derived from the temporal variability measured at a single point by employing the Taylor hypothesis, i.e. when the mean flow transports turbulent eddies past a sensor, angular frequency ω and wavenumber k are related through the mean wind speed u by k = ω/u (Reference TrittonTritton, 1988, p. 308). Assuming that the Taylor hypothesis holds, and that the mean wind at the surface during our measurements was also ∼0.5ms⁻¹, then we estimate a horizontal length scale of ∼0.5 m for pressure variations at a frequency of 1 Hz. Alternating high- and low-pressure perturbations on the scale of decimetres could cause horizontal pressure gradients and shallow horizontal flow in the firn (Fig. 3) on the scale of the sensor depth (0.10 m) instead of vertical pressure gradients and deep penetration as predicted by the one-dimensional treatment. These considerations suggest the need for a three-dimensional theory of wind pumping and motivate the present study.

Fig. 2. Power spectra of air-pressure dato from Agassiz Ice Cap, Arctic Canada. Pressure sensors were placed at the snow surface and at a depth of 0.10 m in firn; output from the two sensors was simultaneously recorded. The wind speed during these obser-vations was 5.1 ms⁻¹ at 0.60 m height. Best-fitting straight lines have been drawn and the pressure spectrum at 0.10m has been predicted using the one-dimensional theory of wind pumping. The discrepancy between the predicted spectrum and the observed spectrum at 0.10 m indicates shortcomings of the one-dimensional theory.

Fig. 3. Penetration of a surface-pressure disturbance into a permeable half-space. If the pressure distribution on the z = 0 plane is spatially variable, as assumed in the three-dimensional theory, then air-flow paths have a lateral as well as a vertical component.

Theory

Penetration of surface-pressure fluctuations into permeable firn

Following the development of Reference Clarke, Fisher and WaddingtonClarke and others (1987), the penetration of a surface-pressure fluctuation p into a uniform permeable half-space is given by

(1)

where µ. is the viscosity of air, n is the porosity of the permeable medium, k is the permeability and p0 is the ambient pressure. We assume that the total pressure is Ρ = po + p where p0 is constant and the fluctuating component p has a mean value of zero. The coordinate system is chosen so that z represents depth, measured positive downward, and the z = 0 plane is located at the boundary between firn and the atmosphere (Fig. 3). Equation (1) is parabolic and its complete solution requires two boundary conditions and an initial condition. As boundary conditions, we assume that p(x, y, 0, t) is known and p(x, y, ∞, t) vanishes. The influence of any initial condition diminishes with time and for our analysis can be neglected.

The Fourier transform of p(x, y, z, t) with respect to x, y and t is written as

(2)

where k_x and k_y are the wavenumbers respectively associated with the x and y spatial coordinates and ω is angular frequency. The corresponding inverse transform is

(3)

Using the derivative properties of Fourier transforms (e.g. Reference BracewellBracewell, 1978, p. 117), Equation (1) can be written

(4)

A general solution of Equation (4) is

(5)

where A(k_x, k_y, ω) and B(k_T, k_y ω) are integration constants. We define

(6)

(7)

and sgn(ω) denotes the algebraic sign of ω.

To obtain β, it is necessary to find the square root of the complex function Υ defined in Equation (6). The standard approach is to decompose Υ into real and imaginary parts, Υ _r = k ² _x + k ² _y and Υ_i = 2sgn(ω) α², and express Υ as the product of modulus and phase terms

(8)

where

(9)

The square root of Equation (8) gives

(10)

From the trigonometric identities,

and

and Equation (9), it follows that

(11)

(12)

(13)

and hence that the real and imaginary parts of β are respectively

(14)

(15)

As previously stated, the boundary conditions that we apply to Equation (5) are that p(x, y, 0, t) is specified and that p(x, y, ∞, t) = 0; equivalent boundary conditions, expressed in terms of Fourier-transformed variables, are that

is specified and that

. These boundary conditions determine the integration constants A(k_x, k_y, ω) and B(k_x, k_y ω). Without loss of generality, we select the positive values of β _r and β _i and, upon applying the boundary conditions, find that Α(k_x, k_y, ω) = 0 and

Thus, the transformed pressure solution is

(16)

Thermal effects of frictional dissipation

Equation (16) shows that rapid fluctuations in surface pressure drive rapid fluctuations in the firn. Any single fluctuation releases a negligible amount of frictional energy; however, the cumulative effect could be important. We examine this by following a stochastic approach that leads to a simple formulation in terms of autocorrelation and power-spectrum functions. Frictional heating is modelled by a distributed heat source a(x, y, z, t) and we assume that firn has constant thermal conductivity K, density ρs and specific heat capacity c. The diffusion of temperature fluctuations is therefore described by the thermal diffusion equation

(17)

If the boundary surface z = 0 is kept at a constant temperature and if the pressure fluctuations causing the frictional heating are assumed to be spatially homogeneous random processes in the x and y variables and stationary in t, then the dependence on x, y and t can be eliminated by spatial and temporal averaging. We shall use ovcrlining to denote averaging over x, y and t variables so that

(18)

The ideas of continuity and differentiability can be extended to random processes (e.g. Reference MiddletonMiddleton, 1960, p. 68; Reference ParzenParzen, 1962, p. 83). In practice, the conditions for differentiability do not place onerous restrictions on T(x, y, z, t) and it is readily shown that

, so that

(19)

The frictional heating term has the magnitude

(20)

where q_i is the air-flux vector and φ is fluid potential. For small pressure fluctuations φ = φ ₀ + p where φ ₀ is large and constant. From Darcy’s law,

so that Equation (20) can be written as

(21)

The time and space average of the frictional heating is therefore

(22)

We define the autocorrelation function of p(x, y, z, t) with respect to x, y and t variables as

(23)

The Wiener-Khintchinc theorem (e.g. Reference MiddletonMiddleton, 1960, p. 141–43) establishes that autocorrelation and power-spectral functions are a Fourier transform pair; thus, the power spectrum of p(x, y, z, t) (holding z untransformed) is given by

(24)

and the inverse transform relation is

(25)

The autocorrelations

and

of the spatial derivatives of p(x, y, z, t) are related to S[p](k_x, k_y, k_z, t) as follows:

(26)

(27)

(28)

Expressed in terms of autocorrelation functions (26)–(28), the averaged frictional heating can be written as

(29)

From the above expressions, it follows that

(30)

Thus the time- and space-averaged temperature equation can be written as

(31)

Equation (31) is a differential equation in z and its solution, subject to the boundary conditions T (0)=0 and

is

(32)

using Equations (7), (14) and (15).

Influence of spectral shape on energy penetration

Insight into the behavior of the temperature solution in Equation (32) can be obtained by substituting assorted simple expressions for S^[p](k_x, k_y, 0, ω) into Equation (32) and performing the necessary integrations. Here the aim is to examine spectral forms that result in tractable integrations when substituted into Equation (32). The simplest assumption is that

(33)

This form corresponds to applying a sinusoidal pressure forcing at the z = 0 boundary. The angular frequency of the fluctuation is ±ω₀ and, at any given time, pressure is constant on the z = 0 plane. The resulting temperature disturbance is

(34)

i.e. that for one-dimensional wind pumping (Reference Clarke, Fisher and WaddingtonClarke and others, 1987).

A more versatile assumption is that the spectral form is described by

(35)

where σ_x and σ_y are constants that determine spectral breadth. A special case of Equation (35) is to take σ_x = σ_y = σ _x so that the pressure fluctuations are spatially isotropic. For this case the spectral function has radial symmetry and it is convenient to define a radial wavenumber k ² _r = k ² _x . + k ². The spectral function can then be written as

(36)

Applying Equation (25), the autocorrelation functions corresponding to Equations (35) and (36) are

(37)

and

(38)

where r ² ₀ = x ² ₀ + y ² ₀. The definitions l_x = 1/σ_x, l_y = σ_y and l ² _x = l ² _x + l ² _y clarify the point that inverse spectral breadths such as 1/σ_x are related to spatial coherence lengths such as l_x . Exploiting the radial symmetry, the solution form in Equation (32) can be written as

(39)

with

(40)

(41)

Performing the ω integration gives

(42)

where α₀ = α(ω₀).

Discussion

The previous section sets out the three-dimensional theory of wind pumping without offering much explanatory detail. Here, we summarize the main results and demonstrate how they can be applied. We shall assume that random pressure fluctuations at the z = 0 boundary between firn and air are described by the power spectrum S ^[p] (k_x, k_y , 0, ω). The penetration of pressure fluctuations into the firn is governed by Equation (16) and, when applied to the power spectrum, gives

(43)

where β _r is evaluated using Equations (7) and (14). Equation (43) can be regarded as an expression that defines a power-spectral transfer function

(44)

between the spectrum of surface-pressure fluctuations and that of fluctuations propagated through firn to a depth z. The main point of interest in Equation (44) is that depth attenuation of the pressure fluctuations is influenced by both the temporal and spatial frequency content of the pressure signal. To illustrate the use of Equation (44), we make the additional simplifying assumption that pressure fluctuations are spatially isotropic so that k ² _x + k ² _y = k _r where k _r is a radial wavenumber. Next, we consider the attenuation of a single-frequency component wu and introduce dimensionless variables z* = z/z ₀ and k* = k _r/k ₀ where z ₀ = l/2α₀, α₀

(45)

which applies to the penetration of a single-frequency component if statistical isotropy in the x and y space coordinates is assumed.

Figure 4 presents a set of graphs of the function ǀH(k*,z*)ǀ² as given in expression (45). The abscissa in Figure 4 is dimensionless wavenumber fc* and the ordinate is the power-transfer function ǀH(k*,z*)ǀ². Labels on individual curves represent values of dimensionless depth z*. We now discuss how this information can be used. Suppose that one wishes to consider the penetration of 0.1 Hz pressure fluctuations into firn. From knowledge of the ambient pressure po, air viscosity μ and firn properties n and κ, the constant α ₀ and the characteristic depth z ₀ can be calculated. If, for example, α ₀ = 0.1 m⁻¹ then z ₀ = 5 m and k ₀ = 0.1 m⁻¹. For the chosen frequency, the attenuation of spatial fluctuations having wavenumbers that are small compared to k ₀ (k* << 1 in Figure 4) will be ruled by the physics of one-dimensional wind pumping, i.e.

(46)

Fig. 4. Transfer function (as a function of dimensionless wavenumber) for the penetration of a surface-pressure disturbance to different depths in the firn. The labels represent the dimensionless depth z* corresponding to each curve. Surface pressure is assumed to vary sinusoidally with time but to have a full spectrum of spatial frequency components so that the forcing at z = 0 is not equivalent to a plane wave of pressure. Note that large values of spatial frequency (large k*) attenuate more rapidly with depth than small values (k* << 1).

whereas spatial fluctuations having wavenumbers large compared to k ₀ (k* >> 1 in Figure 4) will be more strongly attenuated than predicted by the one-dimensional theory. Thus Figure 4 illustrates an obvious, though central, point of this paper: the one-dimensional theory overestimates the penetration depth of real pressure fluctuations, as was also noted by Reference ColbeckColbeck (1989).

Now wc consider the attenuation ǀH(k*,z*)ǀ² predicted by the three-dimensional theory for the pressure signals in Figure 2. Is it compatible with the data? If we assume that the attenuation in Figure 2 is due entirely to three-dimensional effects, then at any frequency we can calculate the horizontal wavenumber k _r and check whether the resulting wavelength 2π/k _r is reasonable. The above estimates α ₀ = 0.1 m⁻¹ and z ₀ = 5.0 m are typical of cold firn such as at Agassiz Ice Cap. With these values, our sensor depth of 0.10 m corresponds to z* = 0.02. Putting the observed attenuation of ǀHǀ² = 0.017 at 1 Hz into Equation (45), wc find a wavenumber of k* = 410 or k _r = 41 m⁻¹. The corresponding wavelength is λ = 0.15 m. The large k* demonstrates that turbulence at 1 Hz is deep into the region where penetration is controlled by three-dimensional geometry.

Eliminating k* from Equation (45) using the Taylor hypothesis (frequency ω and wavenumber k are related through the mean wind speed u by k = ω/u (Reference TrittonTritton, 1988, p. 308)) yields an exponential relation between ǀH(z*ǀ and frequency, rather than the power law implied by Figure 2. This suggests that either the Taylor hypothesis is a coarse approximation at best, or our sensors introduce localized pressure disturbances that render the spatial components of the pressure field non-stationary.

Nonetheless, because these calculated wavelengths are of the same order as those estimated in the Introduction using the observed wind speed, results from Reference BergenBergen (1980), and the Taylor hypothesis, we conclude that three-dimensional effects can account for a large part of the attenuation in Figure 2.

The thermal effects of pressure fluctuations are a second concern of this paper. Expression (32) gives the main result of the three-dimensional theory. As yet, we do not know the actual form of the pressure-fluctuation spectrum S ^[p](k_x, k_y , 0, ω), so there is little profit in examining the general result. We therefore turn attention to Equation (42) which applies to spatially isotropic surface-pressure fluctuations at a single frequency ω ₀. To simplify further discussion, it is helpful to transform Equation (42) to dimensionless variables. We define a dimensionless temperature T* = Τ/T0 where T₀ = _KS ₀/2μΚ and employ the same definitions of z ₀, α ₀ and k ₀ as used previously. Lastly, we introduce a dimensionless parameter A = α ₀/σ _r that characterizes the comparative magnitudes of the horizontal length scale for atmospheric pressure fluctuations and the vertical length scale for pressure penetration. From Equation (38) and the accompanying discussion, it is evident that l ₀ = 1/σ _rindicates the correlation length of atmospheric turbulence cells. Thus the Λ parameter might also be defined as A = l ₀/2z ₀, a ratio of horizontal and vertical length scales. Applying the foregoing scaling to Equation (42) yields the expression

(47)

Fig. 5. Plot of dimensionless temperature Τ* against dimensionless depth z* for a range of values of the parameter A. The dotted curve represents the solution for one-dimensional wind pumping. Note that for large values of A the three-dimensional wind-pumping effect approaches that for the one-dimensional case. The parameter A characterizes the comparative magnitudes of the horizontal length scale for atmospheric, pressure fluctuations and the vertical length scale for pressure pene-tration. The dotted lines indicate particular values z* = 1 and T* = 0.632 which define a boundary-layer thickness for the one-dimensional theory. The depth at which T* = 0.632 decreases for decreasing values of A.

Figure 5 shows the results of evaluating Equation (47). The abscissa is dimensionless depth z* and the ordinate is dimensionless temperature T*; A values are given for each curve and the limiting case of one-dimensional wind pumping is indicated by dashes. We begin by discussing the graph for one-dimensional wind pumping. In terms of the three-dimensional theory, this situation applies when the correlation length of surface-pressure fluctuations becomes infinitely large. As shown in Figure 5, the thermal influence of wind pumping increases with depth. For z* = 1, the dimensionless temperature has increased to a value of T* = 1 − exp(−l) = 0.632, i.e. 63.2% of the full temperature offset T* = 1. (Dotted lines in Figure 5 indicate z* = 1 and T* = 0.632.) We can think of the condition T* = 0.632 as defining the thickness of a thermal boundary layer within the firn. For A = 2.0, the three-dimensional theory predicts a temperature disturbance similar to that for the one-dimensional case, but the depth to the 63.2% temperature influence is reduced. For smaller values of A, the boundary-layer thickness is further reduced. In summary, three-dimensional effects become significant when the correlation length l ₀ is small relative to the characteristic depth z ₀, and the effect is to reduce the thickness of the thermal boundary layer.

From Figure 5, it is apparent that the actual magnitude of the temperature offset is always T* = 1, irrespective of the boundary-layer thickness. In Nature, though not in our theory, the boundary-layer thickness is likely to be important in deciding whether or not wind pumping can significantly disturb the thermal regime of firn at windy sites. If the boundary layer is thin, then during calm periods the firn can re-establish thermal equilibrium with overlying air and there would be no long-term consequence of wind pumping. We believe that boundary-layer thicknesses of the order of ∼5 m are required for significant retention of thermal disturbances. Thus, the three-dimensional theory, which predicts a reduced boundary-layer thickness, also implies that the one-dimensional theory overestimates the thermal consequences of wind pumping.

Concluding Remarks

Determination of palcotcmpcratures from the isotopic record in polar ice cores, for example those taken from Agassiz Ice Cap, requires knowledge of the relationship between 10 m firn temperature and regional mean air temperature. The δ¹⁸O record yields relative changes in cloud temperature at the time of ice-crystal formation (Reference Dansgaard, Johnsen, Clausen and GundestrupDansgaard and others, 1973). Regional palcotemper-aturcs can then be determined from isotope records if snow always accumulates at a uniform rate throughout the year, and if the modern regional average air temperature and average isotopic composition of snowfall can both be measured. Absolute rather than relative values of palcotcmpcraturc are necessary to estimate temperature-dependent processes in the past, such as melt rates on ice caps and sea ice, or range limits of plant assemblages. If, as at Agassiz Ice Cap, the 10 m firn temperature is not a good proxy measurement for mean regional air temperature, then the palcotcmpcratures derived from δ¹⁸O will be uncertain by an un-known constant comparable to the microclimate noise in the 10 m temperatures. Our present study leaves unexplained the cause of anomalous temperatures at Agassiz Ice Cap and raises doubts about the importance of wind pumping as a mechanism for energy transfer to firn.

Glacier ice provides the only continuous samples of the paleoatmosphere (e.g. Reference Lorius, Jouzel, Merlivat, Barkov, Korotkevich and KotylakovLorius and others, 1985). Because ice-core records have been affected by a number of physical processes during deposition and during the long process of compaction into glacier ice, time series of chemical or physical properties measured in ice cores cannot always be simply related to palcoatmospheric conditions. Before we can understand the full climatic significance of ice-core records, we must first understand the physical processes controlling incorporation and possibly alteration of the measured chemical and physical properties with the glacier firn and ice. Wind pumping might prove to be an important process for redistributing atmospheric constituents, stable isotopes and impurities.

Acknowledgements

This research was supported by grants from the Natural Sciences and Engineering Research Council of Canada and by grant DPP-861965 from the U.S. National Science Foundation. Logistical support from the Polar Continental Shelf Project enabled the field measurements leading to Figure 2. We thank J. Cunningham for measurement and processing of the pressure data of Figure 2, and K. Hutter for suggesting improvements to the presentation.