Effects of firn ventilation on isotopic exchange

Abstract

A new model of isotopic diffusion in the upper few meters of firn tracks the isotopic composition of both the ice matrix and the pore-space vapor through time in two dimensions. Stable isotopes in the vapor phase move through the firn by diffusion along concentration gradients and by advection. Wind-driven ventilation carries atmospheric water vapor into the firn, where it mixes with existing pore-space vapor. Unlike previous models, our model allows disequilibrium between pore-space vapor and the surrounding snow grains. We also calculate the isotopic effects of ventilation-driven sublimation and condensation in the firn. Model predictions of isotopic diffusion in firn compare favorably with existing diffusion models. Model results quantify what other investigators have suggested: isotopic change in the upper few meters is more rapid than can be explained by the Whillans and Grootes (1985) model; isotopic equilibration with atmospheric vapor is an important component of post-depositional isotopic change; and ventilation enhances isotopic exchange by creating regions of relatively rapid sublimation and condensation in the firn.

1. Introduction

Stable-isotope ratios recorded in polar snow can be related to the condensation temperature of the original snow grains (Dansgaard, 1964; Jouzel and others, 1997). Variations in temperature are recorded in successive layers of snow as variations in the stable-isotopic ratios of the snow. Through time, isotopic diffusion modifies the isotope profile in snow and destroys high-frequency information. Recent work has focused on reconstructing the original isotopic ratios in snow from the diffused records collected in field samples (e.g. Cuffey and Steig, 1998; Bolzan and Pohjola, 2000). Reconstructing the original records requires an understanding of the processes involved in isotopic diffusion.

Prior work (e.g. Whillans and Grootes, 1985) identified diffusion through the vapor phase as the most important process in redistributing isotopes deep in the firn. In order to predict isotopic changes in snow, prior studies have generally assumed that the pore-space vapor and the surrounding snow grains are in isotopic equilibrium. This assumption allows for an analytical model of isotopic diffusion in the firn, originally presented by Johnsen (1977) and simplified by Whillans and Grootes (1985). Other studies (e.g. Cuffey and Steig, 1998; Johnsen and others, 2000) have modified this model, but have not substantially changed the theory of Whillans and Grootes (1985). Results of these models show broad agreement with observations of isotopic profiles in firn, but do not fully illustrate the underlying physical processes.

Rather than parameterizing isotopic changes in the firn, the model presented here takes a step towards a physically based, time-dependent model of stable-isotopic change in firn. The model tracks the evolution of the isotopic composition of both the pore-space vapor and the firn in a vertical two-dimensional (2-D) cross-section. Isotopes are redistributed in the firn by four processes:

1. (1) Forced airflow (ventilation or wind pumping; Clarke and others, 1987) advects atmospheric vapor (with a specified isotopic composition) into the firn.

2. (2) Isotopes diffuse along isotopic concentration gradients in the pore-space vapor.

3. (3) Vapor and solid phases evolve toward isotopic equilibrium.

4. (4) Sublimation and condensation in the firn alter bulk isotope ratios.

Throughout this work, “condensation” refers to the vapor-to- solid phase change, and “sublimation” refers to the reverse, i.e. solid-to-vapor phase change.

A key innovation in our new model is our explicit treatment of isotopic exchange between phases. Previous authors assumed that pore-space vapor and surrounding snow grains were in isotopic equilibrium; in contrast, our new model allows for isotopic equilibration between the phases over time, without assuming that equilibrium is reached.

2. Post-Depositional Change

2.1. Model overview

Changes in stable-isotope ratios are tracked through time in a 2-D slice of firn. Unlike earlier models of stable-isotope change in firn, which are analytical (e.g. Whillans and Grootes, 1985), ours is an iterative numerical model.

Our isotope model requires estimates of firn temperature and of vapor density in the pore spaces. Both are calculated using the model of Neumann (2003), which is similar to the model of Albert (2002). Isotopic ratios of both the ice grains and the pore-space vapor are calculated on the same control-volume grid (Patankar, 1980) as the temperature and vapor-density fields. Throughout this paper, x and z are the horizontal and vertical coordinates.

Isotopic composition is measured on the δ scale (e.g. Criss, 1999, p.31). The isotopic composition δ_i of the ice matrix is modified by sublimation and condensation in the firn, and by isotopic exchange between the pore-space vapor and surrounding ice grains. The model assumes rapid isotopic homogenization within each ice grain, but neglects long-range isotopic diffusion through the bulk ice matrix, since diffusion through the solid is several orders of magnitude slower than diffusion through the firn. It also neglects diffusion along quasi-liquid surface layers, since the thickness of this layer is very small (~ a monolayer) at temperatures typical on the Antarctic plateau (T < − 30°C).

Three different mechanisms change the isotopic composition δ_v of the vapor:

1. (a) pore-space advective diffusion, which includes processes (1) and (2) in section 1 above, i.e. advection of pore-space vapor and diffusion along isotopic gradients in the vapor,

2. (b) isotopic equilibration between the vapor and the surrounding ice matrix, process (3) above, and

3. (c) phase changes, process (4) above, both sublimation and condensation in the firn.

Before describing the algorithm used to calculate the evolution of δ_v, we outline the procedure used to account for each of these distinct physical processes.

2.2. Pore-space advective diffusion

Atmospheric vapor is advected through the firn by subsurface airflow, which is calculated following Colbeck (1989), as described by Neumann (2003). The magnitude of subsurface airflow depends on the firn microstructure (porosity and permeability) and the often poorly known surface microtopography (Waddington and others, 1996). Airflow through the firn not only carries atmospheric vapor into the firn, but also provides a mechanism to transport pore- space vapor rapidly in the upper 2 m of the firn. In addition, isotopes in the pore-space vapor diffuse along isotopic gradients in the vapor.

With these assumptions, the change in the isotopic composition of the vapor is given by the conservation equation (Appendix A):

(1)

where is the isotopic composition of the vapor under process (a) (pore-space advective diffusion), D _s is the isotopic diffusivity of the pore-space vapor through snow and is the 2-D subsurface air flux (m³ m^-2 s^-1

D _s is based on the diffusivity of water vapor in air; two factors are used to relate D _s to the free-air diffusivity. First, the porosity is used to correct for the fact that in a given volume, only a fraction of the space is occupied by water vapor, the rest being occupied by solid ice. As in other models of isotopic exchange (e.g. Johnsen and others, 2000), we use the total porosity; we do not distinguish between the open and closed porosity discussed by Schwander (1989). The second factor accounts for the blocking effect of ice grains and the path length of channels in the firn. The models of Jean-Baptiste and others (1998) and Johnsen and others (2000) use the tortuosity factor, while the model of Schwander (1989) uses the porosity. Each of these models ultimately relies on the firn density to derive either the porosity or the tortuosity. We use the parameterization of Johnsen and others (2000), which agrees with the revised data of Schwander and others (1988) (personal communication from S. Johnsen, 2003). We do not distinguish between the slightly different diffu- sivities of the heavy and light isotopes of either oxygen or hydrogen. We also neglect the increase of isotopic diffusivity D _s in snow with decreasing air pressure (Cuffey and Steig, 1998), although it could easily be included.

2.3. Isotopic equilibration

Our model allows for isotopic exchange between the pore- space vapor and the surrounding ice matrix. This contrasts with other studies (e.g. Whillans and Grootes, 1985; Johnsen and others, 2000) which assumed that the two phases are in isotopic equilibrium. Following Ingraham and Criss (1998), we use an equilibration function (Equation (2)) in which the condensed phase is an essentially infinite reservoir of the minor isotope, i.e.

(2)

where k is the equilibration rate constant, t is time, the superscript 0 refers to an initial value at t = 0, and the superscript b refers to isotopic composition as a result of pro-cess (b) (isotopic equilibration). The isotopic composition δ_eq of the equilibrium vapor is given by:

(3)

where δ_i is the isotopic composition of the surrounding ice grains and a is the well-known temperature-sensitive fractionation coefficient (Majoube, 1971; Criss, 1999). The factor of 1000 is a result of the definition of the δ scale (Criss, 1999, p.31).

We use Equation (2) to predict changes in isotopic compositions of the solid and vapor in time-steps of Δt. By comparing the initial and final isotopic composition of the vapor, and using information about grid size, porosity and vapor density (which varies exponentially with temperature), we determine the number of heavy and light isotopes exchanged between phases in a given time-step. This procedure calculates the changes in the isotopic composition of both the pore-space vapor and the surrounding ice matrix. Our approach is described in detail in Appendix B. Johnsen and others (2000) outlined a similar approach.

The equilibration constant k has not yet been measured for isotopic exchange between ice and water vapor at any temperature. Our model assumes that the condensed phase represents an infinite reservoir of the minor isotope. Consequently, we assume that k is related primarily to the time required for the pore-space vapor to equilibrate, and this, in turn, is determined by the water-vapor diffusivity in snow.

We now estimate the magnitude of k in two steps. First, Epstein and Mayeda (1953) measured the isotopic equilibration between CO₂ gas and liquid water at 25 ° C. They found that k for this process was ~2.85 h^-1. Based on the ratio of the diffusivities for CO₂ gas and water vapor at 25°C (Massmann, 1998), we estimate that k for isotopic exchange between water vapor and liquid water may be ~5 h^-1. Second, by assuming a linear relationship between k and saturation vapor pressure (where both saturation vapor pressure and k go to zero at absolute zero), we estimate that k may be approximately 1 h^-1 at −25°C.

In the following calculations, we used a range of k from 0.1 h^-1 to 5 h^-1. This estimate neglects the effects on the equilibration constant k arising from the geometry of the condensed phase (Ingraham and Criss, 1993), from the ambient pressure and from the impurity concentration in the condensed phase (Criss, 1999, p. 145). Laboratory measurements are needed to determine an appropriate range of k for isotopic studies at low temperature.

2.4. Sublimation/condensation

The model also includes the isotopic effects of sublimation and condensation in the firn. The phase-change field is calculated as described in Neumann (2003) using the governing equations of Albert (2002). Sublimation and condensation in the firn result from interaction between firn temperature and the local relative humidity in the firn. If airflow through the firn (ventilation) brings cool and relatively dry pore-space air into a warmer area, that pore- space air becomes sub-saturated, and some surrounding ice sublimates. Conversely, if vapor is advected from a warm area into a cooler area, the pore-space air becomes supersaturated, and some vapor condenses onto the surrounding snow grains.

In the isotope model, vapor derived from sublimating snow grains has the same isotopic composition as the solid. In this way, sublimation changes the isotopic composition of pore-space vapor, but does not change the composition of the remaining snow grains directly. However, condensation has isotopic implications for both the solid and the vapor phases. We assume that condensation in the firn follows a Rayleigh process, and the condensate forms in isotopic equilibrium with the vapor, as in Dansgaard (1961):

(4)

where δ_c is the isotopic composition of the condensate, the superscript 0 denotes an initial value at the beginning of the time-step of Δt, the superscript c refers to the isotopic composition of the vapor as a result of process (c) (sublimation or condensation) and β(Δt)is the fraction of the vapor remaining at the end of the time-step Δt. δ_i is updated using a mass-weighted average of δ_c and the isotopic composition δ_i of the pre-existing ice in a unit volume:

(5)

where the superscript 0 denotes an initial value, M _i is the initial mass of ice in a unit volume and M _c is the mass of condensed vapor in a unit volume. The masses of both the ice matrix and vapor are updated accordingly. This model assumes fast diffusion in the ice grains (e.g. Whillans and Grootes, 1985), such that grains have uniform δ_i. This is in contrast to the model of Rempel and Wettlaufer (2003), which tracks changes in the radial distribution of stable isotopes in individual snow grains. M _c is a function of the condensation rate and the time-step Δt. During summer (T ~ −20°C), the time-scale to radially equilibrate grains of radius 500 µm is ~20 days (Whillans and Grootes, 1985); in winter, the time-scale for equilibration is longer (~100 days). Consequently, the assumption of homogeneity is likely to be valid only during summer.

Under these assumptions, the isotopic composition of the remaining vapor is found from (Dansgaard, 1961):

(6)

If pore-space vapor is in isotopic equilibrium with the surrounding snow grains, then small amounts of condensation will not change the isotopic composition of the solid, because the initial condensate will have the same isotopic composition as the surrounding snow grains.

The isotopic composition of the vapor as a result of process (c) (sublimation and condensation) can be summarized as:

(7)

where the superscript 0 again denotes an initial value, the superscript 1 denotes a final value at the end of the time interval Δt, ρ_V is the vapor density in the pore space and s is the sublimation (s > 0) or condensation (s < 0) rate, calculated as in Albert (2002), integrated over the time-step Δt. Neumann (2003) showed that ventilation-driven condensation and sublimation causes only slow growth of snow grains, under polar conditions. Therefore, it is unlikely that grain growth due to condensation will cause large radial isotopic gradients in snow grains. Consequently, we do not include intra-granular diffusion in our model.

2.5. Solution procedure

Equation (1) is solved using the two-dimensional (2-D) control-volume method of Patankar (1980), on the same 2-D grid that we used to calculate the temperature and vapor-density fields (Neumann, 2003). Multiple processes influence both δ_V and δ_i in a given time-step; however, over small time-steps Δt, we can treat them as independent. We solve Equations (1), (2) and (7) using the same initial condition for δ_V and δ_i in each equation. We then add the changes due to each of the three processes, to find the net changes and to determine final values of δ_V and Si at the end of the time- step. This approach is analogous to an explicit scheme (Press and others, 1992), in which changes during each time-step are calculated based on conditions that prevail at the start of the time-step.

This isotopic model requires several boundary conditions in addition to those needed for the heat- and vapor- transport models. We assume that there is no isotopic exchange across the left, right or bottom boundaries of the solution domain, and that the atmospheric vapor has known isotopic composition and vapor density; these may be time- dependent.

The numerical model is run with 1s time-steps. Initial transients in δ_v die out after approximately 5.5 hours. After approximately 5.5 hours, isotopic rates of change are constant in time in both the vapor and the solid. The model is used primarily to calculate the isotopic composition of the solid as a function of time, although other ancillary parameters, such as δ_v; δ_c and the total mass exchanged between phases are calculated as well. In this paper, we present model results for changes in δ¹⁸O. The model could also be used to calculate changes in δD by changing the appropriate parameters (fractionation coefficient α, equilibration constant k and diffusivity D _s)

3. Results

3.1. Effective isotopic diffusivity

In order to maintain analytical simplicity, other models (e.g. Cuffey and Steig, 1998) generally use an effective isotopic diffusivity D _f to account for all of the physical processes, including isotopic diffusion and isotopic exchange between phases, that are included explicitly in our new model. This parameter D _f describes, in an approximate way, how a stable-isotope profile in firn will smooth out over time as a result of these processes.

We can address several questions by comparing our full- physics model to diffusion-only models. First, what values of D _f are most appropriate for diffusion models under various conditions of snow temperature, snow density and ventilation? Second, can we constrain the fractionation rate factor k?

In order to compare the results of our numerical full- physics model with analytical “diffusive” models, we use a simple diffusion-only equation:

(8)

where D _f is the effective isotopic diffusivity. If we assume that there is zero net isotope flux into or out of the ends of the profile, then this equation has a well-known solution:

(9)

where A _n are the Fourier coefficients for the initial profile. We use an initial sinusoidal δ¹⁸O(z, 0) pattern in the firn with a 10% amplitude, mean δ¹⁸O = −40% and an annual wavelength of 0.25 m. This corresponds roughly to conditions on the Antarctic plateau.

By using Equation (9) over the same total elapsed time t, and by varying D _f until the diffusion-only model results and the full-physics model results match as closely as possible, we determine an effective isotopic diffusivity for the full-physics model.

We ran several simplified simulations to calculate stable- isotopic changes in the firn as a result of transport through the vapor phase and isotopic equilibration between solid and vapor, in the absence of convection and temperature gradients. In these runs, firn temperature was uniform at −30°C, airflow velocities in the firn were set to zero, and there was no exchange as a result of vapor-density gradients (because there is no net condensation or sublimation in isothermal snow). The model ran for 2 × 10⁴ s, with a time- step of 1s and vertical resolution of 0.01 m.

Using the best-guess value for the isotopic equilibration coefficient (k = 3 h^-1; section 2.3) leads to an effective isotopic diffusivity D _f in the firn of 1.1 × 10^-5 m² a^-1. This estimate is smaller than the value used by Johnsen and others (2000) (D _f = 5 × 10^-5 m² a^-1) for equivalent temperature (T = − 30°C) and snow density (ρ = 600kgm^-3).

As discussed above, there is considerable uncertainty in the value of k. We repeated this calculation, using a range of values for k in the model, and found the corresponding “best” D _f to represent the model results. Figure 1 shows the dependence of D _f on k. The solid line indicates the most likely range of k for isotopic equilibration between water vapor and snow with density ρ = 600 kg m^-3.

Fig. 1. Effective isotopic diffusivity D _f (m² a^-1) as afunction of isotopic equilibration coefficient k (h^-1). Firn temperature is uniform (= −30^oC) p _snow = 600 kgm ³, and there is no air motion, and no condensation or sublimation. We expect 0.1 ≤ k ≤ 5 h^-1, as indicated by the solid line.

Diffusivity D _s of pore-space vapor is primarily a function of snow density p which is used to determine both the porosity and the tortuosity, as in Johnsen and others (2000). In these simplified model runs, we used a number of different values of p for a given temperature (T = −30°C) and equilibration rate constant (k = 3 h^-1), and found that the resulting effective diffusivity D _f (expressed in m² a^-1) can be expressed as a linear function of p (expressed in kg m ³) (r² = 0.98):

(10)

3.2. Isotopic evolution in summer

We use the results of Neumann (2003) as estimates of the temperature, vapor-density and mass-exchange fields in the firn during summer. The Neumann (2003) model solves the coupled equations for heat and water-vapor motion in a 2-D cross-section of the firn, using the governing equations of Albert (2002). In the model results presented below, we use the approximate steady-state temperature and mass- exchange pattern in the firn corresponding to mid-summer conditions with moderate winds, as given by Neumann (2003).

Airflow in the firn results from steady wind (5 ms ¹ at 10 m height) over a sinusoidal surface topography typical of the Antarctic interior (wavelength of surface topography = 3 m, bump height = 0.1 m), as in Waddington and others (1996). Airflow velocity vectors, calculated as in Colbeck (1989), are shown in Figure 2. Maximum air velocity in the upper few cm of firn is ~0.25 cm s, and airflow is relatively unimportant below 21 m depth.

Fig. 2. Airflow in the firn (ventilation) results from steady wind (5 m s^-1 at 10 m) over a sinusoidal surface topography (wavelength λ = 3 m, amplitude h = 0:1. Following Colbeck (1989), air-pressure variations are applied on a plane that approximates the actual sinusoidal snow surface. Cunningham and Waddington (1993) showed that this is a reasonable approximation if h/λ«1.

Firn temperature from Neumann (2003) is shown in Figure 3. In regions of air inflow (strongest inflow on left and right margins, as shown in Fig. 2), the near-surface firn is close to isothermal as a result of the vertical advection of relatively warm air from the atmosphere (T _atm ~ −20°C). In regions of air outflow (in the center, as shown in Fig. 2), vertical advection carries relatively cold air from the firn interior towards the surface. Horizontal advection between these two regions helps to create the curved pattern of isotherms in the upper few meters. Below 2 m depth, advection is less important, and the temperature field is dominated by conduction, resulting in horizontal isotherms.

Fig. 3. Near-surface firn temperature (in ºC) with subsurface airflow as in Figure 2 and environmental conditions typical of Antarctic summer. Near-surface, near-isothermal zones exist where relatively warm air from the atmosphere (T _atm ~−20ºC) enters the firn. Cold air carried up toward the surface in regions of outflow produces large lateral temperature gradients. Below 2 m, temperature is dominated by conduction.

The phase-change field from Neumann (2003) is shown in Figure 4. Advection of air through temperature gradients in the firn produces regions of supersaturated and sub-saturated air in the pore spaces. In regions of air inflow, warm saturated air enters the firn. This air is supersaturated relative to the cooler firn at depth, resulting in regions of condensation (negative mass exchange centered at ^0.5 m depth). In regions of air outflow, relatively cool and dry air from the firn interior is advected into progressively warmer firn near the surface. This air is sub-saturated, resulting in sublimation of some of the surrounding snow. Lateral advection between regions of inflow and outflow, coupled with the non-horizontal isotherms, curves the locus of maximum sublimation toward the surface.

Fig. 4. Mass-exchange rate calculated using environmental conditions typical of Antarctic summer with subsurface airflow as in Figure 2. Positive values indicate sublimation; negative values indicate condensation. Air flowing into the firn reaches saturation vapor pressure in the upper 0.05m, regardless of the relative humidity of the atmosphere (Neumann, 2003). Advection of relatively warm moist air from the atmosphere into progressively colder and drier firn results in condensation. Sublimation occurs in outflow regions, as cold dry air from the firn is advected into progressively warmer near-surface regions.

The temperature and phase-change fields are treated as constant with respect to time. These fields provide an idealized “snapshot” of the firn during mid-summer. In order to generate an equivalent snapshot of isotopic changes in the firn during mid-summer, we ran the isotope model for 20 000 s (time-step = 1 s, vertical resolution = 0.01−5 m, horizontal resolution = 0.06−4 m).

We assume that the isotopic composition δ_atm of the near-surface atmospheric vapor is constant and uniform, and is in isotopic equilibrium (Equation (3)) with the surface (i.e. summer) snow (Massey, 1995). Air flowing into the firn reaches saturation vapor pressure in the upper few cm (Neumann, 2003); consequently, for tracking mass and isotopic changes at depth, we neglect this near-surface mass exchange and assume that air flowing into the firn is at saturation vapor density. In fact, this near-surface sublimation (or condensation) may alter the isotopic compositions of both the upper cm of snow and the vapor entering the firn. Consequently, available data on atmospheric relative humidity should be used when applying this model to a specific site.

Airflow velocities in the firn are negligible below ~2 m with the representative microtopography in our model. We expect that below the upper ~2 m, isotopic change in the firn can be well described by an analytic diffusion model, such as Equation (8).

Predictions of the isotopic change of the firn in regions of maximum air inflow are shown by the solid line in Figure 5. Below 2 m depth, the model predicts that changes in δ_i follow a standard diffusion model; both summer and winter snow layers are modified equally, and the mean change is zero. The lower solid line with dots represents a fit of a diffusion-only model (Equation (9)) to the numerical model results. At ~2.5 m depth, the effective diffusivity D _f of the numerical model is 2.4 × 10^-5 m² a^-1. D _f decreases with increasing depth, as shown from 2 to 3 m in Figure 5, because the reduced firn temperature lowers the vapor pressure. At lower vapor pressures, there are fewer water molecules in the vapor phase, so the rate of isotopic exchange between phases decreases.

Fig. 5. Numerical-model predictions of isotopic change of firn over 5.5 hours in regions of air inflow with environmental conditions typical of Antarctic summer. Atmospheric vapor is in equilibrium with surface (summer) snow (δ_atm ~ −43%). Equilibration with atmospheric vapor causes asymmetric changes in near-surface layers (winter snow modified more than summer layers). Below ~2 m, only diffusion in the vapor phase and equilibration between solid and vapor drive isotopic changes. Dashed line and solid line with dots show best fit of diffusive model (Equation (9)) to model with full physics, to determine effective dffusivity Df of full-physics model.

It is apparent that other processes are active in the upper 1.5m of the firn in Figure 5. In the upper 30 cm, summer layers are not noticeably modified, while winter layers are strongly modified, producing an asymmetrical pattern. Investigating the isotopic composition δ_v of the pore-space vapor (Fig. 6) reveals that, in regions of air inflow, δ _v is strongly influenced by δ_atm We assume that δ_atm is in equilibrium with the uppermost summer layers, and by extension, far out of equilibrium with the winter layers. As this vapor is advected into the firn, isotopic equilibration causes a large change in the winter layers, but relatively small change in the summer layers.

Fig. 6. Numerical-model prediction of pore-space vapor (δ_v, solid line with dots), firn (δ_i, solid line) and equilibrium vapor (from Equation (3), dashed line). In the upper 1m, isotopic composition of pore-space vapor is strongly influenced by atmospheric vapor (δsub>atm ~ -43‰). Below 2 m, δ_v is controlled by a balance between diffusion in the vapor phase and equilibration between solid and vapor.This model predicts that pore-space vapor is rarely in isotopic equilibrium with the surrounding snow, in contrast to the equilibrium assumption of Whillans and Grootes (1985).Vapor disequilibrium is a result of rapid diffusion in the vapor phase, and relatively slow isotopic equilibration between solid and vapor (Equation (2)).The equilibration constant k must be ~100 h^-1 in order to generate pore-space vapor in equilibrium with the solid, a value two orders of magnitude higher than our expected range of k (Fig. 1).

Condensation makes only relatively minor contributions to changes in the isotopic composition δ _v of the vapor in regions of air inflow. In the firn, the maximum condensation rate in regions of air inflow (Fig. 4) is found at ~ 0.5 m depth. Pore-space vapor, above ~ 1 m, has advected recently from the atmosphere, and its isotopic ratio is close to the ratio of atmospheric vapor. Isotopic equilibration between this vapor and surrounding snow (in particular, winter snow layers) lowers δ _v significantly faster than the decrease in δ _v due to condensation. By 1 m depth, the pore- space vapor attains an isotopic value that is in equilibrium with a bulk mixture of summer and winter layers, so equilibration of advected moisture has little impact on δ _v. However, condensation preferentially removes heavy isotopes from vapor and decreases δ_v (Equation (2)) regardless of the isotopic composition of the firn, so condensation is a leading cause of isotopic change in the vapor at this depth. The pore-space vapor at 1 m depth is slightly closer to equilibrium with the winter layers, resulting in condensate which is isotopically more similar to the winter layers, leading to greater changes in the summer layers. In addition, the remaining pore-space vapor is now closer to equilibrium with the winter layers than with the summer layers, so isotopic equilibration causes larger changes in the summer layers. However, these two effects are subtle, and are much smaller than the changes in the upper 0.75 m due to advection of atmospheric vapor into the firn. Below ~2 m, condensation rates are very small, and the pore-space vapor becomes more like the vapor in diffusion-dominated systems.

If we represent all the physical processes in the upper firn through a single effective diffusivity D _f, equilibration between atmospheric vapor and the near-surface firn in the upper 1.5 m in regions of inflowing air leads to relatively large values of D _f. The upper dashed line in Figure 5 represents the best fit of the model with only “diffusion” (Equation (9)) to our model with full physics. At 0.5 m depth, the numerical model has an effective diffusivity of 1.4 ×10^-4 m²a^-1; this is 75% higher than D _f at an equivalent temperature in a model in which molecular diffusion is the only process.

The solid line in Figure 7 shows the model prediction of isotopic change in the firn in regions of maximum air outflow. Below ~2 m, these areas are well modelled by a standard diffusion equation; both summer and winter layers are modified equally, and the mean change is 0. The lower solid line with dots in Figure 7 represents the fit of an analytical model (Equation (9)) to the numerical results. At ~2.75 m depth, the effective diffusivity Df of the numerical model is 2.3 × 10^-5 m² a^-1 By ~3m depth, D _f is equivalent between regions of air inflow and outflow. This supports the idea that ventilation effects are not important below ~2.5 m, given the airflow field in Figure 2 resulting from our choice of microtopography.

Fig. 7. Numerical-model predictions of isotopic change over 5.5 hours in regions of air outflow. Asymmetric change in near-surface layers (winter modified slightly more than summer layers) due indirectly to sublimation. Sublimation makes δ_v heavier and closer to isotopic equilibrium with summer layers; consequently, isotopic equilibration changes winter layers more rapidly. Below ~2 m, isotopic changes are due to diffusion in the vapor phase and equilibration between solid and vapor. Dashed line and solid line with dots show best fit of diffusive model (Equation (9)) to model with full physics, to determine effective diffusivity D _f of full-physics model.

In the upper 1 m, the winter layers are modified slightly more than summer layers. The reason for the asymmetry can be found by investigating the pore-space vapor (Fig. 8). Above 1 m, δ_v becomes progressively heavier as a result of the sublimation of surrounding snow grains (Fig. 4). Consequently, δ _v more closely resembles vapor in equilibrium with summer layers. Since δ _v is farther out of equilibrium with winter layers than with summer layers, winter layers are modified slightly more through isotopic equilibration. This effect intensifies towards the surface, where sublimation rates are highest. The upper dashed line in Figure 7 shows the fit of an analytical model to the numerical model results, and reveals that at ~0.5m depth, D _f is 6.7 × 10^-5 m² a^-1. Thus, ventilation increases D _f by ~20% compared to modelled D _f without ventilation effects at an equivalent temperature.

Fig. 8. Numerical-model prediction of pore-space vapor (δ_v, solid line with dots), firn (δ_i, solid line) and equilibrium vapor (from Equation (3), dashed line) in regions of air outflow. Increase in mean δ_v between 3 and 2 m is due to the temperature sensitivity of the fractionation coefficient α (Majoube, 1971). Above 2 m, sublimation of surrounding snow grains (Fig. 4) slightly accelerates the increase in δ_v, because we assume that snow grains sublimate without fractionation. Increased airflow near the surface smooths out variations in δ_v there.

The increase in the mean δ_v value between 3 and 2 m in Figure 8 is due to the temperature sensitivity of the fractionation coefficient α (Majoube, 1971). This effect is present throughout the entire depth, as a result of the vertical temperature gradient (Fig. 3). Above 1 m, the influence of sublimating snow grains causes greater changes in δ _v. The increase in mean δ_v as a result of the temperature dependence of α is also present in Figure 6 below 1.5 m, although condensation (which decreases δ_v) obscures this dependence above 1.5 m.

3.3. Isotopic evolution in winter

As a result of the non-linearity of water-vapor pressure over ice (Colbeck, 1990), near-surface ρ _v is lower in winter by a factor of 25 compared with summer vapor density, and mass-exchange rates are as much as two orders of magnitude smaller in winter than in summer (Neumann, 2003). Model results predict that D _f is an order of magnitude smaller in winter than in summer. Consequently, this analysis neglects the influence of winter ventilation on stable-isotope exchange.

4. Discussion

These results suggest that summer is the most important season for post-depositional isotopic exchange in the firn, as suggested previously (e.g. Waddington and others, 2002). Since temperatures are higher in summer, isotopic exchange is more rapid. However, if the firn models of Albert (2002) and Neumann (2003) generate realistic estimates of subsurface mass exchange and temperature in the firn, then isotopic exchange during winter and other seasons should not be ignored. In particular, condensation in the firn results in potentially important isotopic exchange during any season. The condensation rates reported here result from moderate winds (5 ms^-1 at 10 m height) and modest surface microtopography. Stronger winds will increase condensation rates in the firn, and increase the depth to which advection is important (Neumann, 2003).

This model currently uses a single temperature and vapor-pressure field, because recalculating these fields through time is computationally intensive. Model results suggest that mass exchange (condensation and sublimation) plays a minor role, and temperature changes will primarily affect the vapor density in the firn. Changes in vapor density change the number of water molecules in the vapor phase available for isotopic exchange, and so influence the effective isotopic diffusivity directly. The model results presented here also show that the advection of atmospheric vapor is an important component, which suggests that changes in ventilation strength through time (including wind speed) and the isotopic composition of the atmosphere (δ_atm) could be important. However, it is apparent that the isotopic response time of the firn is slow, and that the firn will not be sensitive to high-frequency changes in ventilation strength or temperature. Consequently, we favor using a monthly averaged temperature field and wind speed in our model.

If we assume that the above model is representative of isotopic exchange during the summer months, we can estimate what the total change might be during a summer season. As noted above, winter layers are modified more rapidly than summer layers during the summer. At ^0.5 m depth in regions of air inflow, this model predicts that winter layers change by 10^-3‰ over 0.23 days. This rate of change would result in a change of ~0.25% over 60 days. In outflow regions, the rate of change is about 50% smaller than in inflow regions, leading to a change of ~0.1 % during summer. This assumes that the isotopic composition of the atmospheric vapor is constant. In fact, the isotopic composition of the atmosphere is strongly influenced by weather systems, and can change by ~10‰ in a day (Grootes and Stuiver, 1997). However, it is apparent that ice grains in the firn respond slowly to isotopic changes in atmospheric vapor. Consequently, we suggest that firn will respond to the average isotopic composition of the vapor over a period of months, and will not be strongly influenced by high-frequency changes in the isotopic composition of atmospheric vapor.

For these environmental conditions, our model predicts that full isotopic equilibration between near-surface winter layers and atmospheric vapor will not be reached. This is supported by the model of Waddington and others (2002) which was based on the characteristic times of the isotopic exchange processes involved. For the environmental conditions used here (T ~ − 30°C; 10 m wind speed = 5ms^-1; accumulation rate = 0.25 m a^-1) the model of Waddington and others (2002) predicts that isotopic equilibrium between the firn and the atmospheric vapor will not be reached at any depth, but that some change toward equilibrium may occur. As noted by Waddington and others (2002), firn will never actually reach isotopic equilibrium with atmospheric vapor, since δ_atm changes more rapidly than the firn can equilibrate, effectively changing the “target” towards which the firn is equilibrating.

The rate-limiting step for isotopic exchange in the firn is isotopic equilibration between the pore-space vapor and the surrounding ice grains. Our model predicts that S _v is far from equilibrium with the surrounding snow grains (Figs 6 and 8). Larger values of k in Equation (2) promote more vigorous exchange between solid and vapor, leading to larger values of D _f. However, full isotopic equilibration between firn and the atmosphere would require unrealistically large values of k. Other models (e.g. Rempel and Wettlaufer, 2003) assume isotopic equilibrium at the boundary between the solid and fluid phases (an essentially infinite k between the pore-space vapor and the outermost shell of the ice grain, but with Equation (2) modified to account for diffusion in the solid as the rate-limiting step). The model of Rempel and Wettlaufer (2003) uses a single representative value for the ice-grain radius. By considering (among other processes) equilibration with firn, which is an aggregate of grains with a variety of sizes and shapes, the model presented here obscures the importance of isotopic diffusion within individual grains in favor of a bulk equilibration factor k. We have estimated a range for k based on available data, and we anticipate that our model will be useful in firn, which generally has a wide range of grain-sizes.

Our model currently provides a snapshot of rates of change in the firn at a single time. If the necessary boundary conditions are known (e.g. surface temperature, relative humidity of the atmosphere, and δ_atm) the model could be used in transient applications. Currently, model application is limited by this lack of data.

As a result of this snapshot approach, the model cannot accommodate the addition of new snow at the surface. As other investigators have shown (e.g. McConnell and others, 1997; Van der Veen and others, 1999), snow accumulation in Antarctica is episodic, and it can occur at any time of the year. As snow accumulates at the surface, near-surface layers are advected down through zones of sublimation or condensation and equilibration in the upper 2 m. Consequently, accumulation of snow at the surface changes the spatial relationship between the isotopic profile in the snow (which is advected downwards relative to the snow surface) and regions of rapid phase change (which depend on surface topography and wind speed, but remain at a fixed depth relative to the snow surface).

In areas with high accumulation rate (e.g. much of Greenland), snow is advected rapidly through the ventilated zone. Consequently, we anticipate that the rapid isotopic exchange in the upper meter described here will not be important in such areas. In such areas, a computationally simpler model such as Johnsen and others (2000) can be used. The effects of firn ventilation on isotopic composition described here become more important as accumulation rate decreases, as described by Waddington and others (2002). In particular, we expect to observe the effects described here in areas where the accumulation rate is ⋍ 5 cma^-1 ice equivalent or less, such that snow remains in the ventilated zone for many years.

Modelling the ongoing compaction and densification in the firn will also be an important component of an improved model of post-depositional isotopic change in the firn. Our model uses a uniform density profile in the firn, and a uniform wavelength of δ(z). As noted by several investigators (e.g. Johnsen, 1977; Rempel and Wettlaufer, 2003), strain decreases the wavelength of isotopic signals through time. Decreasing the wavelength of S increases the spatial gradient of the pore-space vapor, leading to greater diffusive fluxes. In addition, metamorphic processes (such as grain boundary sliding and grain growth) act to increase firn density ρ_snow through time. Increases in ρ _snow cause decreases in D_f through more efficient blocking of vapor flow (Cuffey and Steig, 1998).

Several model parameters are poorly known. As discussed above, there is large uncertainty in k. In addition, uncertainty in D _s could arise from errors in D _a, which is ultimately based on work presented by Geiger and Poirer (1973) and used by several studies (e.g. Whillans and Grootes, 1985; Cuffey and Steig, 1998; Johnsen and others, 2000). Our model also may underestimate vapor flux in the firn, since the influence of water-vapor flow due to radius-of- curvature grain growth (e.g. Colbeck, 1980) is not considered.

Although isotopic diffusion in the firn is a slow process, the cumulative effects can be large. In low-accumulation- rate areas, snow remains in the near-surface ventilated zone where isotopic exchange is relatively rapid for many years. It seems likely that winter snow will be modified more than summer snow, leading to an apparent truncation of the winter portion of the annual S cycle in snow. If this is correct, it may not be possible to distinguish this type of postdepositional isotope change from the wind-scour model of isotopic change proposed by Fisher and others (1983). Fisher and others noted that an isotope record from a windy divide site at Agassiz Ice Cap, Ellesmere Island, Canada, was apparently missing much of the winter precipitation, based on comparisons with the isotope record from the adjacent flank site, which is much less windy. In the isotope model presented here, stronger winds will lead to larger D _f by increasing the flux of water vapor from the atmosphere into the snow. Equilibration with atmospheric vapor will tend to reduce the amplitude of winter layers, which is also the effect of the wind-scour model of Fisher and others (1983). However, it is unlikely that post-depositional processes are responsible for a significant fraction of the isotopic change at Agassiz Ice Cap, due to the relatively large annual accumulation rate there (~10−20 cm a^-1 ice equivalent; Fisher and others, 1983). Unless accounted for, either wind scour or the post-depositional asymmetric isotopic change suggested here could lead to erroneous reconstruction of the magnitude of seasonal precipitation or temperature variations, particularly at sites with low accumulation rates.

5. Conclusions

This work represents a step towards a process-based understanding of post-depositional stable-isotope change in the firn. The model results confirm what other investigators have suggested: that isotopic smoothing in the upper few meters is more rapid than can be explained by the Whillans and Grootes (1985) model, particularly in low-accumulation-rate areas; that isotopic equilibration with atmospheric vapor is an important component of postdepositional isotopic change (Waddington and others, 2002); and that firn ventilation can enhance isotopic exchange (Neumann, 2001).

If ventilation effects are neglected (no airflow, sublimation or condensation), the effective diffusivity D _f of this model broadly agrees with other models (e.g. Johnsen and others, 2000). This model suggests that isotopic equilibrium between the pore-space vapor and surrounding snow grains is rarely achieved, in contrast to the assumption of equilibrium in other parameter-based models such as Cuffey and Steig (1998). The success of those models may be due in part to tuneable parameters such as critical density (Cuffey and Steig, 1998) or the application of those models to high-accumulation-rate areas (Johnsen and others, 2000) where ventilation effects should be less important. Our model lacks some mechanisms for vapor exchange (such as grain growth) that would increase Df and improve agreement with Johnsen and others (2000) (see section 3.1). The disagreement could also be due in part to the large uncertainty in k (Equation (2)). It is unlikely that k is large enough to allow equilibrium between phases. Another approach (taken by Rempel and Wettlaufer, 2003) is to assume equilibrium between the ice surface and the surrounding pore- space vapor, and to treat diffusion within the ice grains as the rate-limiting step.

Future laboratory work will focus on reducing the uncertainty in k, while future fieldwork should be undertaken to measure δ_atm. In future model developments, we plan to incorporate snow accumulation, grain metamorphism and strain. In addition, isotopic inhomogeneities within individual snow grains may be important, and their influence on isotopic diffusion in the firn should be examined further.