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. Reference Whillans and GrootesWhillans 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 Reference NeumannNeumann (2003), which is similar to the model of Reference AlbertAlbert (2002). Isotopic ratios of both the ice grains and the pore-space vapor are calculated on the same control-volume grid (Reference PatankarPatankar, 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. Reference CrissCriss, 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 Reference ColbeckColbeck (1989), as described by Reference NeumannNeumann (2003). The magnitude of subsurface airflow depends on the firn microstructure (porosity and permeability) and the often poorly known surface microtopography (Reference Waddington, Cunningham, Harder, Wolff and BalesWaddington 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. Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others, 2000), we use the total porosity; we do not distinguish between the open and closed porosity discussed by Reference Schwander, Oeschger and LangwaySchwander (1989). The second factor accounts for the blocking effect of ice grains and the path length of channels in the firn. The models of Reference Jean-Baptiste, Jouzel, Stievenard and CiaisJean-Baptiste and others (1998) and Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) use the tortuosity factor, while the model of Reference Schwander, Oeschger and LangwaySchwander (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 Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000), which agrees with the revised data of Reference Schwander, Stauffer and SiggSchwander 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 (Reference Cuffey and SteigCuffey 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. Reference Whillans and GrootesWhillans and Grootes, 1985; Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others, 2000) which assumed that the two phases are in isotopic equilibrium. Following Reference Ingraham and CrissIngraham 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 (Reference MajoubeMajoube, 1971; Reference CrissCriss, 1999). The factor of 1000 is a result of the definition of the δ scale (Reference CrissCriss, 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. Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen 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, Reference Epstein and MayedaEpstein 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 (Reference MassmanMassmann, 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 (Reference Ingraham and CrissIngraham and Criss, 1993), from the ambient pressure and from the impurity concentration in the condensed phase (Reference CrissCriss, 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 Reference NeumannNeumann (2003) using the governing equations of Reference AlbertAlbert (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 Reference DansgaardDansgaard (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. Reference Whillans and GrootesWhillans and Grootes, 1985), such that grains have uniform δ_i. This is in contrast to the model of Reference Rempel and WettlauferRempel 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 (Reference Whillans and GrootesWhillans 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 (Reference DansgaardDansgaard, 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 Reference AlbertAlbert (2002), integrated over the time-step Δt. Reference NeumannNeumann (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 Reference PatankarPatankar (1980), on the same 2-D grid that we used to calculate the temperature and vapor-density fields (Reference NeumannNeumann, 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 (Reference Press, Teukolsky, Vetterling and FlanneryPress 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.1. Effective isotopic diffusivity

In order to maintain analytical simplicity, other models (e.g. Reference Cuffey and SteigCuffey 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 Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen 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 Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen 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 Reference NeumannNeumann (2003) as estimates of the temperature, vapor-density and mass-exchange fields in the firn during summer. The Reference NeumannNeumann (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 Reference AlbertAlbert (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 Reference NeumannNeumann (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 Reference Waddington, Cunningham, Harder, Wolff and BalesWaddington and others (1996). Airflow velocity vectors, calculated as in Reference ColbeckColbeck (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 Reference ColbeckColbeck (1989), air-pressure variations are applied on a plane that approximates the actual sinusoidal snow surface. Reference Cunningham and WaddingtonCunningham and Waddington (1993) showed that this is a reasonable approximation if h/λ«1.

Firn temperature from Reference NeumannNeumann (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 Reference NeumannNeumann (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 (Reference NeumannNeumann, 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 (Reference MasseyMassey, 1995). Air flowing into the firn reaches saturation vapor pressure in the upper few cm (Reference NeumannNeumann, 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 α (Reference MajoubeMajoube, 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 α (Reference MajoubeMajoube, 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 (Reference ColbeckColbeck, 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 (Reference NeumannNeumann, 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.