List of Symbols
1. Introduction
The stable-isotope ratios δ 18 O and δD measured in ice from polar cores provide detailed records of past temperature (Reference DansgaardDansgaard,1964; Reference JohnsenJohnsen and others, 1997; Reference PetitPetit and others, 1999). The isotope records themselves are smoothed by diffusion at a rate that is enhanced by connected vapor pathways in the firn (Reference JohnsenJohnsen, 1977; Reference Whillans and GrootesWhillans and Grootes, 1985). Several inversion methods have been developed to account for this post-depositional alteration and recover the original isotopic signals (Reference JohnsenJohnsen,1977; Reference Cuffey and SteigCuffey and Steig,1998; Reference Bolzan and PohjolaBolzan and Pohjola, 2000). Seasonal variations are of particular interest, both for their importance to our understanding of climate dynamics and for their use in absolute dating. Beneath the level of pore close-off, post-depositional alteration is slowed considerably. Nevertheless, diffusion continues to smooth the isotope records, leading to the loss of seasonal and other high-frequency information. In fact, Reference JohnsenJohnsen and others (1997) found that, in the deep Greenland ice sampled by the Greenland Icecore Project (GRIP) core, the rate of smoothing is considerably faster than would be predicted by diffusion through solid ice alone. Intriguingly, this apparent excess diffusion appears to apply only within the Holocene ice at GRIP; the rate of alteration within ice from the Wisconsin glacial period is consistent with the normal solid-state diffusivity. Reference NyeNye (1998) argued that the presence of liquid veins at the boundaries of ice grains (Reference Nye and FrankNye and Frank, 1973; Reference MaderMader, 1992a, Reference Maderb; Reference Nye, Maeno and HondohNye,1992) might explain the excess diffusion. However, the analysis of Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) suggests that the required vein dimensions are unrealistically large. Here, we model the diffusion of stable isotopes in polycrystalline ice and show that the predictions of Reference NyeNye (1998) and those of Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) represent two end-members of a more general range of behavior.We identify the physical characteristics that determine which of these asymptotic regimes more closely resembles the prevailing conditions and quantify the role of liquid in the smoothing of isotopic signals. A deeper understanding of the processes by which the isotope records evolve will help to increase the accuracy and resolution of paleotemperature reconstructions.
In polycrystalline ice that is below the normal bulk-melting temperature of 0°C, veins of pre-melted liquid trace the triple junctions between three grains and join at nodes where four ice grains meet (Reference Nye and FrankNye and Frank, 1973). The veins and nodes constitute a connected liquid network that is prevented from freezing at sub-zero temperatures by two physical effects. First, the surface energy of the solid–liquid interface favors the presence of liquid as the curvature of the vein walls increases; this is the Gibbs–Thomson effect. Second, soluble impurities that are rejected from the ice lattice are concentrated in the liquid and reduce its chemical potential so that it remains in equilibrium with the solid at colder temperatures. Intermolecular interactions between adjacent ice grains or between the ice and a foreign substance such as a gas bubble or a dust particle are a third cause of pre-melting that can also contribute to the total quantity of liquid present (Reference Ohtomo and WakahamaOhtomo and Wakahama, 1983; Reference Dash, Fu and WettlauferDash and others, 1995; Reference WettlauferWettlaufer, 1999; Reference Wettlaufer and DashWettlaufer and Dash, 2000). The liquid network facilitates the transport of soluble impurities and assists the purification of temperate ice (Reference Glen, Homer and ParenGlen and others, 1977; Reference MaderMader,1992b). In polar ice, the diffusion of soluble impurities through the liquid network can influence the compositional records recovered from ice cores (Reference Rempel, Waddington, Wettlaufer and WorsterRempel and others, 2001, Reference Rempel, Wettlaufer and Waddington2002; Reference Rempel and WettlauferRempel and Wettlaufer, 2003) and enhance the potential for chemical reactions to take place (Reference FisherFisher, 1987). Here we examine how pre-melted liquid affects the stable-isotope records beneath the firn.
Reference NyeNye (1998) proposed his model, referred to here as the Nye model, immediately after Johnsen discussed his observations at a NATO Advanced Study Institute in 1997 (Reference Johnsen, Clausen, Jouzel, Schwander, Sveinbjörnsdóttir, White, Wettlaufer, Dash and UntersteinerJohnsen and others,1999). In the Nye model, diffusion through the vein liquid is treated as instantaneous and the rate of diffusive smoothing is determined by the exchange of isotopes between the solid grains and the liquid veins. Adopting parameter values that represent conditions typical of the GRIP Holocene ice, the model predictions reproduce the order-of-magnitude increase in apparent diffusivity that Reference JohnsenJohnsen and others (1997) inferred from the data. Existing inversion algorithms, which might be used to recover the original isotope record, cannot be applied directly to the Nye model. Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) addressed this issue by developing an alternative formulation for vein-mediated isotopic diffusion, referred to here as the Johnsen model, which produces predictions that are at odds with the Nye-model results. In particular, using parameter values that correspond to those adopted by Nye, theJohnsen model predicts a more modest three-fold increase to the apparent diffusivity. TheJohnsen model is based on earlier treatments of isotopic diffusion in the firn (Reference JohnsenJohnsen, 1977; Reference Whillans and GrootesWhillans and Grootes, 1985), with the liquid network playing the role that the vapor phase takes higher up in supplying fluid pathways for rapid diffusion.The ice grains and the vein liquid at each depth are assumed to be in local isotopic equilibrium and the rate of diffusive smoothing is controlled by diffusion through the fluid phase. It is worth noting here that the rate of signal alteration in the Nye model is independent of the liquid diffusivity, while the Johnsen-model predictions are insensitive to the diffusivity of the solid. The model we present here simulates the entire process including both the exchange of isotopes between the phases and the rapid transport through the fluid. We begin by presenting a basic model that illuminates the essential physics using the geometry considered by Reference NyeNye (1998). We then generalize the model to simulate the conditions sampled by ice cores, and examine the evolution of an arbitrary initial signal that is subject to the effects of vertical strain. Informed by comparisons between the model predictions and ice-core observations, we discuss possible causes of the apparent excess diffusion and suggest directions for further research.
2. The Model
To provide physical insight, we seek an analytical description of isotopic diffusion in the model system shown in Figure 1. A cylindrical annulus of ice with outer radius b, chosen to correspond to the grain radius (typically b = O (10−3) m), envelops a fluid cylinder of radius a that represents a vein (typically a = O (10−6) m). The radially averaged volumetric concentratio
of the minor isotope (18O or D) displays an initial periodic signal with wavelength λ in the axial direction. This averaged signal can be thought of as representing the measured seasonal variations in the isotope content of ice-unity. To predict the evolution of the isotope profile, we proceed by examining the conservation conditions in each phase.
Fig. 1. The geometry used to model the vein-mediated diffusion of stable isotopes through polycrystalline ice (see text for details).
In the solid outer cylinder, we consider an arbitrary control volume Ωs with surface S s, shown schematically in Figure 1. Diffusion across the surface of the control volume changes the number of molecules of the minor species contained within it, so the evolution of the number density of the minor isotopic species N′ is described by
where D is the molecular diffusivity and the subscript s refers to the solid. The divergence theorem is applied to the right side and we note that the integration volume is arbitrary so that
where the radial r and axial z coordinates are indicated in Figure 1. A no-flux condition holds along the outer boundary so that
The condition along the inner boundary where r = a is found by matching with the isotope concentration in the fluid phase, to be examined next.
The vein is sufficiently narrow that the concentration in the fluid depends only on z. Accordingly, we consider the conditions for isotope conservation in a disk-shaped control volume Ωv of height Δ and radius a + Δ that spans the vein cross-section, as shown in Figure 1. We confine our discussion to a purely diffusive parameter regime; scaling arguments that elucidate the conditions under which fluid flow through the vein system is an important transport mechanism are provided in the Appendix. Changes in the isotope content are accomplished by diffusion down the axial fluid-concentration gradient and exchange with the solid ice across the phase boundary at r = a so that
where the subscript v refers to the vein liquid.The final term in Equation (4) accounts for the diffusive exchange between the ice and the vein fluid that is driven by gradients in the isotopic content of the solid
at the radial boundary of Ωv. The divergence theorem is applied to the first surface integral and each term is integrated as Δ tends to zero so that, to leading order in Δ, Equation (4) becomes
Isotopic equilibrium between the fluid and the solid along their common boundary implies that
where N represents the concentration of the major species (16O or H) and the fractionation coefficient α is close to unity We substitute Equation (6) into Equation (5) and use Equation (2) to eliminate the time derivative so that
This is the second boundary condition that solutions to Equation (2) must satisfy.
In ice-core studies the isotope concentrations are normally reported as δ values, which are defined in terms of the standard ratio S according to
Noting that D v ≫ D s, α ≈ 1 and the ratio N s/N v is approximately equal to the ratio of the solid and liquid densities ρ s/ρ v ≈ 1, we transform variables and write the governing Equation (2) and its boundary conditions (3) and (7) as
and
The volume fraction of liquid is much less than unity so that b ≫ a and the radially averaged initial condition can be written as
where kz = 2π/λ is the wavenumber of the axial isotopic variations, and δ 0 and δ 1 are the mean value and amplitude of the initial signal. Equations (9–12) constitute a complete description of isotopic diffusion in the model system.
We use separation of variables and write the solution to Equation (9) that satisfies Equations (10–12) as
where the radial wavenumber k r is found as the solution to the following transcendental equation:
The radial dependence of the isotopic composition is expressed in terms of the Bessel functions J and Y as
which averages to
Equation (13) predicts that the presence of veins in polycrystalline ice causes a harmonic signal with wavenumber kz
to be attenuated at a rate equivalent to that in a homogeneous medium with an effective diffusivity that is
times greater than the molecular diffusivity in an ice monocrystal D
s. The enhancement factor f provides a convenient measure of the effect of pre-melted liquid on the diffusion of stable isotopes. Equation (13) can also be written in terms of a diffusion length
that describes the smoothing effects so that
Signals with wavelengths that are comparable to or smaller than the diffusion length are smoothed significantly, whereas those with much longer wavelengths retain their initial form.
2.1. Model predictions
2.1.1. The Nye model
It is useful to compare our predictions with those of the Nye and Johnsen models. The Nye model proceeds from the observation that D v ≫ D s and examines the limiting case in which the fluid diffusivity may be treated as infinite so that the flux boundary condition given in Equation (11) is replaced by δ| r = a = δ 0. The predicted evolution of δ can be expressed in the same manner as the model presented here, except with the radial wavenumber k r in Equations (13) and (15) replaced by k N, which satisfies
which is the natural consequence of D
v ≫ D
s in Equation (14). Significantly, Equation (17) indicates that k
N is insensitive to the temperature-dependent diffusivities D
v and D
s, and the axial wavenumber kz
. The Nye model can also be written in the form of Equation (16), with the diffusion the length modified so that
, where the corresponding enhancement factor is defined as
.
2.1.2. The Johnsen model
The Johnsen model treats the case in which radial variations in isotope concentration are negligible and all transport takes place through the fluid. For the purposes of comparison with our model results, this second assumption is relaxed to allow axial concentration gradients to drive diffusion through the solid as well as the fluid. The appropriate governing equation is readily obtained from considerations similar to those that led to Equation (4), but with the integrating volumes encompassing a thin slice through both cylinders so that
The assumption of equilibrium between the phases in the Johnsen model implies that F
J(r) ≡ 1 and the isotopic profile is written in the form of Equation (16), but with the diffusion. length given by
, wherein the enhancement factor is f
J = 1 + ϕD
v = D
s and ϕ = a
2/b
2 is the porosity. By including the effect of diffusion through the solid phase along the axial concentration gradient, described by the second term on the right of Equation (18), we have ensured that the diffusion length σ
J does depend weakly on the solid-state diffusivity and the enhancement factor approaches unity at zero liquid fraction. The Johnsen model accounts for misorientations between the vein axis and the direction of the prevailing isotopic gradient by replacing the vein diffusivity with D
v/τ
f, where τ
f is the tortuosity factor. For the modeled geometry, τ
f = 1; later, we show the behavior predicted by our model for an isotropic distribution of vein orientations so that D
v is reduced by τ
f = 3.
2.1.3. Model comparison
We adopt parameter values that are representative of the so conditions sampled by ice cores. The outer radius b is chosen to correspond to a typical grain-size so that b = 10−3 m. The inner radius a is chosen to be typical of expected vein sizes so that a = 10−6 m. The average of the measured diffusion parameters reported by Reference RamseierRamseier (1967) is used to calculate the temperature-dependent molecular diffusivity from
where T is the absolute temperature. For the liquid diffusivity, the experimental data reported by Reference Gillen, Douglass and HochGillen and others (1972) for supercooled water at temperatures down to 242 K are fit to a quadratic given by
The temperature in the upper part of the GRIP core is approximately 241 K; using this as a nominal value, the diffusivities are D s ≈ 9.6×10−17 m2 s−1 and D v ≈ 1.8×10−10 m2 s−1.
The annual-layer thickness decreases with depth, and accumulation rates vary widely, so the axial length scale for isotopic variations is highly sample-dependent. Figure 2 shows the predictions of the three models for the enhancement factor f as a function of wavelength. This is the factor by which the apparent diffusivity of the solid–liquid mixture exceeds that of single ice crystals. For isotopic signals with short wavelengths, our model predictions (solid line) follow those of the Nye model (dot-dashed line), which increases from unity by an amount proportional to the square of wavelength. At longer wavelengths, the enhancement tends to the horizontal asymptote predicted by the Johnsen model (dashed line). The dotted line shows the predicted enhancement when the liquid diffusivity is reduced by a tortuosity factor of τ f = 3.
Fig. 2. The enhancement factor f as a function of wavelength λ = 2π/kz for the present model (solid line), the Nye model (dot-dashed line) and the Johnsen model (dashed line). The dotted line shows the predicted enhancement when the liquid diffusivity is reduced by a factor of τf = 3. The temperature T = 241 K, the vein radius a = 10−6 m and the grain-size b = 10−3 m.
A closer examination of the transcendental Equation (14) for k r reveals that the dimensionless parameter
is diagnostic of the regime with which the model predictions are most closely associated. Physically, ϒ can be thought of as a measure of the relative efficiency of isotopic transport through the fluid compared to isotopic exchange between the phases. We expect the radial wavenumber k r to be largely controlled by the grain-size so that k r = O(b −1) and to leading order in a = b Equation (14) can be approximated as
At small wavelengths, ϒ is large (i.e. ϒ ≈ 7×102 for λ = 0.01 m using the nominal parameter values described above) and Equation (22) reduces to Equation (17) derived for the Nye model. At large wavelengths, ϒ is more modest in size (i.e. ϒ ≈ 7 for λ = 0.1 m) andthe behavior predicted by Equation (22) tends closer to the Johnsen-model results.
Further intuition into the physical processes that control the model behavior is gained by examining the radial isotopic profiles, shown in Figure 3. In the Johnsen model (dashed line) equilibrium is assumed so the isotope concentration is equal to the radial average throughout and F J ≡ 1. In the Nye model (dot-dashed line) the diffusivity of the fluid is infinite so the isotope concentration at the vein wall is equal to the vertical average and F N(a) = 0. The region near the vein wall is expanded in the figure inset, which shows that the present model (labeled solid lines) predicts an isotope concentration at the vein wall (vertical dotted line) approaching unity for λ = 0.1 m. When the length scale of axial isotopic variations is relatively large, the concentration gradient that drives diffusion through the fluid is small, and exchange between the solid and liquid easily replenishes the isotopic content of the veins so that it stays near the radially averaged value. At a wavelength of λ = 0.01 m our model profile is more similar to that predicted by the Nye model. The shorter wavelength produces larger axial concentration gradients that deplete the isotope concentration in the veins more rapidly than diffusion through the solid can resupply it; to compensate, the isotope concentration at the vein wall drops so that the axial gradient is reduced and the radial gradient through the solid is increased. Further analysis shows that the enhancement factor can be written as f ≈ 1 + ϕF(a)D v/D s when Equation (22) holds; this form highlights the connection between our model predictions and those of theJohnsen model, for which F(a) = 1.
Fig. 3. Normalized isotopic profiles
as a function of radial position.The inset shows a close-up near the vein wall (vertical dotted line). Model profiles (solid lines) are shown for signal wavelengths of λ = 0.1 m and λ = 0.01 m, labeled in the inset.The Nye (dot-dashed line) and Johnsen (dashed line) model radial profiles are independent of λ. Calculations were performed using the nominal parameter values of T = 241 K, a = 10−6 m and b = 10−3 m.
Figure 4 shows that all three models are sensitive to the grain-size b. For the Johnsen model, a decrease in b implies an increase in the volume fraction of liquid ϕ = a 2/b 2 and the efficiency of axial transport is increased. In the Nye model, a drop in b reduces the distance for transport through the solid so that the enhancement rises. The present model embodies both of these effects; demonstrating the dominance of radial transport for short-wavelength signals (large ϒ) and changes in liquid volume fraction for long-wavelength signals (smaller ϒ). In addition to affecting the porosity and the distance for radial transport, changes in grain-size alter the value of the dimensionless parameter ϒ. Holdingall else constant, ϒ increases with b so we expect smaller grain-sizes to tend increasingly toward the Johnsen-model predictions and larger grain-sizes to promote behavior that approaches closer to the Nye-model results. Since the Nye model and the Johnsen model each neglect part of the transport process, they both represent upper bounds on the rate of isotopic diffusion predicted by our model.
Fig. 4. The enhancement factor as a function of the grain-size b (T = 241 K and a = 10−6 m).Predictionsforthepresent model and the Nye model are shown for wavelengths of λ = 0.1 and 0.01 m.TheJohnsen model is insensitive to wavelength.
At the temperatures present in polar cores, the volume fraction of pre-melted liquid is expected to be dominated by the effect of impurities in lowering the equilibrium melting temperature (Reference Rempel, Wettlaufer and WaddingtonRempel and others, 2002). For typical impurity loadings and grain-sizes of order 10−3 m, this implies that vein radii are of order 10−6 m. Impurity concentrations change over length scales that are comparable to λ and are correlated with changes in grain-size (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997); to satisfy equilibrium requirements (Reference MaderMader, 1992a; Reference Nye, Maeno and HondohNye, 1992) the vein radii are expected to remain relatively uniform over much larger distances. In Figure 5, the enhancement predicted by the Johnsen model is shown to increase at a rate proportional to a 2 since this is the degree by which the liquid fraction ϕ increases with vein radius. The Nye model is only modestly affected by changes to a, since an increase to the vein radius has a proportionately modest effect on the radial distance for isotopic transport and exchange between the phases. The enhancement predicted by the current model is tracked in the behavior of the Johnsen model at small vein radii, which correspond to smaller values of ϒ. Our predicted enhancement approaches the Nye-model result asymptotically at larger vein radii, as shown by the solid curve for λ = 0.01 m.
Fig. 5. The enhancement factor as a function of the vein size a (T = 241 K, b = 10−3 m, λ as labeled).
2.2. A generalized model
The isotope records sampled by ice cores are not simple harmonic signals, but instead are characterized by variations that span a range of wavelengths. The initial isotopic profile at the base of the firn, immediately after pore close-off, can be expressed as a sum of M harmonic components
with initial wavenumbers kzi
0 = 2π/λi
0, wavelengths λi
0, and amplitudes δ
ic and δ
is. In addition to diffusive smoothing, the effects of vertical strain must be accounted for if inversion techniques are to be used to recover the original isotopic profiles. For a uniform vertical strain rate
, the equation that describes the evolution of the isotopic profiles is
In the ice sheets, recrystallization processes act to ensure that the effects of strain do not significantly alter the geometry of grains as a function of depth (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997). We assume that the diffusive exchange of isotopes between the solid and the liquid is sufficiently rapid that horizontal strain has a negligible effect on the radial isotopic profiles. In this limit, the boundary conditions given by Equations (10) and (11) apply to the solution of the generalized model as well. The solution to Equation (24), subject to Equations (10), (11) and (23), is
where the ith axial wavenumber is kzi
= 2π/λi
and the evoltion of its wavelength is described by
Each harmonic component is attenuated to a degree that depends on the diffusion length σ
i for that portion of the signal, where
with the ith enhancement factor defined as
and each k
ri
satisfying Equation (14).The normalized radial profile for each component Fi
(r) is described by Equation (15).
For physical insight, it is useful to compare the model behavior to that predicted in the limiting cases where the enhancement factor is adequately described by either the Nye model or the Johnsen model. Taking the case of a constant vertical strain rate, theJohnsen-model enhancement factor is substituted into Equation (26) to predict a diffusion length from
which is independent of wavelength. At long times, the exponential term tends to zero
and the diffusion length approaches a steady state, with the tendency for diffusion to increase σ exactly balanced by the effects of vertical strain. In the Nye-model regime, the enhancement factor is inversely proportional to
and the diffusion length for a constant vertical strain rate is given by
where the final term is written in terms of the initial wave-number so that the time dependence is made explicit. At long times, the exponential terms decay and the diffusion length approaches the value that would be predicted for a homogeneous medium with diffusivity D s.
Figure 6 shows the evolution of the diffusion length calculated from Equation (26) with a vertical strain rate of
(chosen to correspond to
, where the accumulation rate i
and the ice-sheet thickness is H = 3000 m). For a harmonic signal with an initial wavelength of 0.01 m, the generalized model agrees with the predictions of the Nye model, with the diffusion length increasing rapidly to a maximum of just over 0.5 cm, before dropping slightly at later times. For an initial wavelength of 0.1 m, the predicted diffusion length tracks theJohnsen-model results at first, then decreases as strain reduces the wavelength and puts the system into the Nye-model regime. The tendency for the diffusion length to reach a maximum at intermediate times rather than increase monotonically to a constant asymptote is a significant qualitative difference between our model predictions and the behavior expected for a homogeneous single-phase medium. Vertical strain eventually reduces the wavelength sufficiently that the Nye-model regime is reached, and at later times the diffusion length tends to the asymptote that would be predicted for a homogeneous medium with the diffusivity of solid ice (dotted line). This behavior is consistent with the observation of a large excess diffusion in the Holocene ice at GRIP, whereas the annual signals in the Wisconsin ice are attenuated much more slowly.
Fig. 6. The evolution of the diffusion length for the generalized model (solid lines), the Nye model (dot-dashed lines) and the Johnsen model (dashed lines) at the labeled initial wavelengths. The dotted line shows the diffusion length predicted for a homogeneous medium with diffusivity Ds. All calculations performed using the nominal parameter values of T = 241 K, b = 10−3 m and a = 10−6 m at a constant strain rate of
.
The attenuation of harmonic signals is shown as a function of their initial wavelengths in Figure 7. Results are presented in terms of their amplitudes after a time interval of 103 years, relative to their initial amplitudes. The predictions of the Nye model follow ours for short wavelengths. The Nye-model predictions tend to a horizontal asymptote of approximately 0:38 at larger wavelengths; the Johnsen-model results are more closely identified with those of our model in the large-wavelength regime. The relative amplitudes predicted by our model are always somewhat higher than those predicted by both the Johnsen and Nye models. The dotted line labeled f = 1 shows the attenuation that would be achieved in a homogeneous medium with the solid-state monocrystalline diffusivity D s. The observations from GRIP Holocene ice reported by Reference JohnsenJohnsen and others (1997) are consistent with a diffusivity that is approximately 30 times greater than D s; the dotted line labeled f = 30 shows the signal strength that is calculated using this enhancement factor. Figure 8 compares the predicted signal attenuation with the attenuation that would be expected for a homogeneous medium with the solid-state monocrystalline diffusivity D s, after the labeled time periods. Also shown are predictions in which the veins are assumed to be randomly oriented with respect to the prevailing isotopic gradient, so that the effective liquid diffusivity is D v/τ f, where the tortuosity factor is τ f = 3.
Fig. 7. The predicted amplitudes of harmonic signals after 1000 years, measured relative to their initial amplitudes and plotted as a function of their initial wavelengths. The dotted lines show the amplitudes that would be predicted for constant enhancements of f = 1 and f = 30. All calculations performed using T = 241 K, b = 10−3 m, a = 10−6 m and
.
Fig. 8. The predicted amplitudes of harmonic signals after the labeled time intervals, measured relative to their initial amplitudes and plotted as a function of their initial wavelengths. The solid lines show the amplitudes predicted using the nominal parameter values.The dashed lines show the dicted amplitudes when the liquid diffusivity is reduced by a tortuosity factor of τf = 3. The dotted lines show the amplitudes that would be predicted for a homogeneous medium with the solid-state monocrystalline diffusivity.
3. The Causes of Excess Diffusion
The model we have presented clearly demonstrates the importance of pre-melted liquid in enhancing the rate of isotopic diffusion in polar ice. Its effect on the relative amplitudes of isotopic signals is illustrated in Figure 8 by the differences between the dotted lines and the solid and dashed lines that represent our model predictions for τ f = 1 and τ f = 3. However, Figure 7 shows that the predicted amplitudes for conditions similar to those encountered in the Holocene ice at GRIP are considerably higher than those observed by Reference JohnsenJohnsen and others (1997), which are represented by the dotted line calculated using their inferred enhancement of f = 30 (Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others, 2000). The predicted enhancement can be increased by decreasing the model grain-size (see Fig. 4) or increasing the model vein size (see Fig. 5), but it is difficult to explain such rapid attenuation without invoking a much higher liquid fraction than is commonly believed. We have considered only the contribution from the bulk colligative effects of soluble impurities and have not considered that intermolecular forces can cause liquid films to form on grain boundaries. Although such grain-boundary melting takes place on a smaller length scale than accessed using typical optical methods, the large surface area involved can result in a significant volume of interfacial pre-melted liquid. In addition to promoting enhanced rates of diffusive transport through the veins, a higher liquid fraction would tend to increase the importance of advective fluid transport, as discussed in the Appendix. Other sources of the apparent discrepancy between the model predictions and ice-core observations require further consideration, discussed presently.
3.1. Liquid geometry
The model is formulated for an idealized geometry that does not capture all the subtle features that characterize the pre-melted-liquid network in polycrystalline ice. To test the sensitivity of the model predictions to the grain-scale geometry, we examined an inverted model in which a cylinder of solid ice with radius b is enclosed by a thin film of liquid with thickness d. The inner radius b is identified with the grainsize, and the film thickness d is chosen so that the liquid fraction ϕ ≈ 2d/b matches that of our original model. For signals with longer wavelengths, the predictions of the inverted model and the original model are identical. This reflects the asymptotic behavior of our model at small ϒ, which yields predictions that are in agreement with those of theJohnsen model and depend on the liquid volume-fraction. At shorter wavelengths, in the Nye-model regime, the inverted model and the original model differ slightly in detail, but show the same qualitative behavior. As noted earlier, our predicted enhancement is bounded above by that of both theJohnsen and Nye models. Since the idealized model geometry does not cause a reduction in the Johnsen-model enhancement, which is itself much lower than the enhancement inferred from the GRIP Holocene ice, we conclude that there must be another reason why the model predicts a weaker enhancement than is observed.
3.2. Grain-size variations
The model is formulated for polycrystalline ice with a constant grain-size and vein size. Theoretical considerations (Reference Nye, Maeno and HondohNye, 1992) and laboratory observations (Reference MaderMader, 1992b) indicate that the vein size should be fairly uniform. The same is not true for grain-sizes, which are anticorrelated with the impurity loading and vary over length scales that are commensurate with those of the annual layering (Reference Thorsteinsson, Kipfstuhl and MillerThorsteinsson and others, 1997). The effects of this variation can be explored with a continuum model that includes variations in liquid volume fraction ϕ. Consider the evolution of the average isotope concentratio
in a control volume that is large in comparison to the grain-size b, but small in comparison to λ. The isotope concentration in the liquid is
, where the scaling factor 1 > F
v > 0 is determined by the geometry of the liquid network; in the idealized geometry considered above, F
v = F(a). To leading order in ϕ, the evolution of the isotopic profiles is given by
where the enhancement factor is f = 1+ϕF v D v/(τ f D s), and the tortuosity factor τ f is included explicitly. Variations in grain-size cause variations in the liquid volume fraction ϕ and scaling factor F v, and lead to gradients in the enhancement factor that produce the third term on the right side of Equation (29). This additional term alters the character of the solution to Equation (29) so that it cannot be written in the same form as we have used for our earlier model. Gradients in f cause different portions of a harmonic signal to diffuse at different rates so that the signal shape is altered in a complicated manner. Depending on their spatial dependence, variations in grain-size could either enhance or reduce the amplitudes inferred for the annual isotopic signals. Further efforts to quantify this effect will require a careful comparison of the spatial relationships between grain-size variations and changes in isotope concentration. At the present stage, a parametric study within the general framework described here would extend beyond the level constrained by firm observational evidence of these controlling effects. Nevertheless, in its current form, the model provides useful guidance for data analysis and interpretation.
3.3. Physical properties
One further possibility that must be considered is whether errors in the physical properties that are used by the model, in particular the diffusivities of the solid and the liquid, might be partly responsible for the apparent diffusion excess. The diffusivity D v was measured in a pure, supercooled liquid (Reference Gillen, Douglass and HochGillen and others, 1972), whereas the pre-melted liquid is in equilibrium, mainly due to the presence of impurities. Impurity ions in solution are surrounded by hydration shells in which the mobility of the water molecules is reduced. Hence, the presence of impurities should tend to make the effective liquid diffusivity lower than that found in supercooled liquid. This suggests that the predicted enhancement can be reduced by the presence of impurities, but it is unlikely to be increased.
The diffusivity of the solid D s was measured on single ice crystals (Reference RamseierRamseier, 1967). Solid-state diffusion in ice involves the motion of point defects, primarily vacancies, so one might expect processes that alter the concentration of point defects to also produce changes to D s. We are not aware of any systematic measurements of D s in ice that has been strained, confined at elevated pressures, or doped with the impurity molecules, NH3, F−and Cl−,thatareknowntoformsubstitutional defects in the ice lattice. Evidence to support the view that the diffusivity measured by Reference RamseierRamseier (1967) might be lower than that which applies to the polar ice-core record comes from a consideration of isotopic diffusion in the firn. Models for firn diffusion (Reference JohnsenJohnsen, 1977; Reference Whillans and GrootesWhillans and Grootes, 1985) have been remarkably successful (Reference Cuffey and SteigCuffey and Steig, 1998) despite their assumption of isotopic equilibrium between the solid and vapor phases. As Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others (2000) note, and a simple adaptation of our model shows, the exchange of isotopes by diffusion between the solid and vapor would be too slow to make this assumption valid if the Reference RamseierRamseier (1967) diffusivity D s applies.While it is possible that recrystallization processes enhance the isotopic exchange in the firn and explain the rapid equilibration (Reference Johnsen, Clausen, Cuffey, Hoffmann, Schwander, Creyts and HondohJohnsen and others, 2000), a higher diffusion coefficient for the solid is an alternative explanation. We note that, for all but the signals with the shortest wavelengths, our model predicts that the solid and liquid beneath the firn are nearly in isotopic equilibrium so recrystallization processes cannot increase enhancement above that which we predict.
4. Conclusions
The stable-isotope records that are recovered from ice cores represent invaluable archives of past temperatures. The accuracy of the paleothermometry that derives from such records relies on a proper treatment of processes that affect the isotope profiles during their long residence times in an ice sheet. The model that has been presented here quantifies the role that diffusion through pre-melted liquid plays in enhancing the rate of isotopic smoothing beneath the firn. The process can be subdivided into two steps: (1) exchange of isotopes between the solid grains and the pre-melted liquid, and (2) rapid diffusion through the liquid along the dominant isotope-concentration gradient. Both processes must be accounted for in a complete description of the diffusive smoothing of isotopic signals, and each comes to the fore asthe dominantcontrol onthe rate of diffusive smoothing in a different parameter regime. When grains are larger and wavelengths are shorter, the lateral exchange of isotopes is the slower, rate-limiting process. In this parameter regime, the isotopic content of the solid grains is laterally inhomogeneous, emphasizing the requirement that sample volumes contain many grains in order to provide representative measurements. When the grains are smaller and wavelengths are longer, pre-melted liquid is most effective at enhancing the rate of diffusion so that isotopic exchange is rapid and the speed of axial transport controls the diffusive smoothing. In this parameter regime, the efficiency of axial transport increases in proportion to the volume fraction of liquid that is present. This implies that variations in impurity loading, which are anticorrelated with variations in grain-size, should produce variations in the rate of axial transport and alter the harmonic content of isotopic signals. This behavior may be responsible for the apparent diffusivity inferred from the GRIP Holocene ice being larger than our model predicts. Alternative explanations that we have explored include the possibility that the liquid fraction is larger than is commonly believed, or that, due to defectscale phenomena, the solid-state diffusivity in the ice sheets is larger than the laboratory-measured value.
Acknowledgements
We are extremely grateful to S. Johnsen andJ. F. Nye for their valuable reviews that have substantially shaped the present paper and have provided suggestions for future work.We note in particular that a correspondence with J. F. Nye over a period of several months was responsible for the inclusion of the Appendix. We thank S. Jones for handling this paper as Scientific Editor. The manuscript has further benefited from the comments and criticisms of J.G. Dash. A.W.R. acknowledges useful discussions with T. Neumann and T. Creyts. Our research is supported by the U.S. National Science Foundation (OPP9908945), Yale University and the Leonard X. Bosack and Bette M. Kruger Foundation.
Appendix
Comparing Advective and Diffusive Transport in Veins
This appendix resulted from a series of exchanges with J. F. Nye during the review process. The principal aspect of these exchanges that we develop here deals with quantifying the relative importance of advection and diffusion in the redistribution of stable isotopes through the vein system in ice sheets.
A comparison of terms in the advection–diffusionequation that describes transport through the veins leads us to define the Peclet number
to measure the relative size of these effects. Here, v is the average velocity of the water in the veins, λ/2π is a characteristic length scale for changes in isotopic concentration, and D v/τ f is the effective molecular diffusivity for a tortuosity factor of τ f ≈ 3. The advective flux is negligible in comparison to the diffusive flux when Pe ≪ 1. Here we estimate the range of wavelengths λ < λ c that satisfy this condition.
In principle, Darcy’s law can be used to obtain the transport velocity U ≡ vϕ, but there are no direct measurements of the permeability Π or the liquid pressure gradient in polar glacial ice. Reference Nye and FrankNye and Frank (1973) discuss an idealized vein network in which they express the permeability as Π = ϕ 2(2b)2/χ where the proportionality constant is χ = 2000. Although Reference Raymond and HarrisonRaymond and Harrison (1975) derive a more general model (see also Reference PhillipsPhillips, 1991, p. 29–34) that accounts for vein obstructions in temperate ice, the Nye and Frank model provides us with a useful upper bound on Π for a given liquid fraction ϕ. We assume that the pressure gradient driving flow is produced by the difference in density between the vein liquid ρ v ≈103 kg m−3 and ice ρ s ≈ 9.2×102 kg m−3 so that Darcy’s law gives U = Π/ν (ρ v – ρ s)g, where g ≈ 10 m s−2 is the acceleration of gravity, and the viscosity at −30°C is ν ≈ 9×10−3 Pa s (Reference Cho, Urquidi, Singh and RobinsonCho and others, 1999). We set the Peclet number to unity and find that the characteristic wavelength λ c that divides the two transport regimes is
Realistic parameter values as used throughout this paper are for b = 10−3 m, ϕ = 10−6 and a liquid diffusivity of D v ≈1.8×10−10 m2 s−1, yielding λ c ≈ 2 m. We conclude that the purely diffusive model presented here should be valid for the parameter regime experienced in a wide variety of ice-core settings.
