Articles | Volume 17, issue 7
Research article
26 Jul 2023
Research article |  | 26 Jul 2023

Isotopic diffusion in ice enhanced by vein-water flow

Felix S. L. Ng

Diffusive smoothing of signals on the water stable isotopes (18O and D) in ice sheets fundamentally limits the climatic information retrievable from these ice-core proxies. Past theories explained how, in polycrystalline ice below the firn, fast diffusion in the network of intergranular water veins “short-circuits” the slow diffusion within crystal grains to cause “excess diffusion”, enhancing the rate of signal smoothing above that implied by self-diffusion in ice monocrystals. But the controls of excess diffusion are far from fully understood. Here, modelling shows that water flow in the veins amplifies excess diffusion by altering the three-dimensional field of isotope concentration and isotope transfer between veins and crystals. The rate of signal smoothing depends not only on temperature, the vein and grain sizes, and signal wavelength, but also on vein-water flow velocity, which can increase the rate by 1 to 2 orders of magnitude. This modulation can significantly impact signal smoothing at ice-core sites in Greenland and Antarctica, as shown by simulations for the GRIP (Greenland Ice Core Project) and EPICA (European Project for Ice Coring in Antarctica) Dome C sites, which reveal sensitive modulation of their diffusion-length profiles when vein-flow velocities reach  101–102 m yr−1. Velocities of this magnitude also produce the levels of excess diffusion inferred by previous studies for Holocene ice at GRIP and ice of Marine Isotope Stage 19 at EPICA Dome C. Thus, vein-flow-mediated excess diffusion may help explain the mismatch between modelled and spectrally derived diffusion lengths in other ice cores. We also show that excess diffusion biases the spectral estimation of diffusion lengths from isotopic signals (by making them dependent on signal wavelength) and the reconstruction of surface temperature from diffusion-length profiles (by increasing the ice contribution to diffusion length below the firn). Our findings caution against using the monocrystal isotopic diffusivity to represent the bulk-ice diffusivity. The need to predict the pattern of excess diffusion in ice cores calls for systematic study of isotope records for its occurrence and improved understanding of vein-scale hydrology in ice sheets.

1 Introduction

The water stable isotope records (δ18O, δD) from ice cores are key proxies for reconstructing palaeoclimatic temperature. Isotope diffusion, which occurs rapidly in firn (mainly by vapour diffusion in the pores) and slowly in ice below the firn transition, causes progressive smoothing that reduces the high-frequency content of these records, strongly attenuating the amplitude of short signals and limiting the depth to which annual and seasonal information survives (Johnsen, 1977; Whillans and Grootes, 1985; Cuffey and Steig, 1998; Johnsen et al., 2000). In the ice cores from central Greenland, annual signals often persist for  104 years, to the early Holocene or just into the glacial part of the record (Johnsen, 1977; Johnsen et al., 1997, 2000), whereas in the ice cores from the East Antarctic plateau, they rarely penetrate through the firn into the ice, owing to low accumulation rates (causing short δ cycles, which decay quickly) and substantial noise and intermittency on the signals during deposition (Laepple et al., 2018; Casado et al., 2020).

Post-depositional diffusive smoothing limits the time resolution of real climate signals extractable from different depths of an ice-core isotope record. The smoothing rate needs to be known in several analyses: (i) studies that use “back diffusion” or deconvolution (Johnsen, 1977; Johnsen et al., 2000) to restore the original annual δ cycles at the surface, for inferring detailed climatic variations (e.g. Küttel et al., 2012; Zheng et al., 2018) or aiding the identification of annual layers in ice-core dating (e.g. Hammer et al., 1978; Vinther et al., 2006); (ii) studies that use the diffusion length σ estimated from the frequency spectrum of isotopic signals at different depths (including where annual cycles are no longer visible) to determine past changes in σ at the firn transition and hence reconstruct the surface temperature history by using the temperature dependence of firn isotope diffusion (e.g. Gkinis et al., 2014; Holme et al., 2018; Kahle et al., 2021); and (iii) studies that model the down-core profile of the diffusion length to assess the climatic variability on different timescales on a record (e.g. Jones et al., 2017; Gkinis et al., 2021; Grisart et al., 2022). The diffusion-length theory of Johnsen (1977), which tracks how σ evolves as a result of isotopic diffusion and vertical compression in the ice column, forms the basis of all these studies.

Whereas the smoothing process in firn has been studied sufficiently to yield models for temperature reconstructions (in (ii) above) and models able to reproduce the observed signal decay in firn (e.g. Cuffey and Steig, 1998; Johnsen et al., 2000; Gkinis et al., 2021), the smoothing process in ice remains poorly understood. We extend its theoretical description in this paper. Diffusive smoothing in ice can strongly impact deep isotopic signals given their long residence time; also, the diffusion rate increases in the warmer ice towards the ice-sheet base. Thus, a key concern motivating our work is that an inaccurate model for the signal smoothing rate in ice can bias the aforementioned studies – notably studies of types (ii) and (iii) above, when applied to records far under the firn transition. This problem may extend to reconstructions that use the differential diffusion length between oxygen and deuterium as a temperature proxy (Simonsen et al., 2011; Holme et al., 2018).

From the decay of annual δ18O cycles along the Holocene part of the GRIP (Greenland Ice Core Project) ice core, Johnsen et al. (1997) inferred an isotopic diffusion rate about 10 times faster1 than the self-diffusion rate in ice monocrystals (Ramseier, 1967), measured at the temperature of the GRIP ice under analysis. They suggested the grain interfaces in polycrystalline ice as causing “excess diffusion” – an idea which prompted three mathematical models seeking to explain the phenomenon. Nye (1998) modelled the effect of water veins located at triple junctions of grain boundaries (Nye, 1989; Mader, 1992a, b) and showed how rapid (liquid phase) diffusion in the vein network “short-circuits” slow diffusion in the ice grains to enhance the signal decay rate above that due to solid diffusion. For signals at the decimetre scale, and ice with a mean grain size of several millimetres, his model predicts an enhancement that matches the GRIP observations. Johnsen et al. (2000) considered more generally interstitial water at grain boundaries as well as in the veins, and calculated how much these pathways raise the effective isotope diffusivity of the bulk ice. Their model couples the isotope concentrations in the solid and liquid in a less sophisticated way than Nye's treatment but accounts for the tortuosity of the veins, and they highlight the possibility for the acidity of the ice to affect the amount of interstitial water. Lastly, Rempel and Wettlaufer (2003), building upon Nye's (1998) continuum description of the grain–vein system, showed that the perfect short-circuiting assumed in Nye's model overestimates the excess diffusion: the enhancement is less than the value predicted by Nye since the liquid diffusivity is high but finite. Rempel and Wettlaufer (2003) clarified the two previous models as end-member approximations of the system, and like Nye's result, their solution gives the enhancement as a function of signal wavelength. All three models – of Nye, Johnsen et al., and Rempel and Wettlaufer – predict a higher enhancement for thicker veins, because wider liquid pathways promote short-circuiting.

Despite these models' implication that recrystallisation and impurity processes in polycrystalline ice can alter the amount of excess diffusion to shape the signals on isotope records in complex ways, no diffusion-length-based models or temperature reconstructions have yet incorporated their results into calculations, which often assume the monocrystal diffusivity of Ramseier (1967) below the firn transition. Nor has there been progress in unravelling the controls and mechanisms of excess diffusion – by theory or experiment – for 2 decades. Yet, the potential occurrence of excess diffusion continues to concern ice-core studies. When analysing the WAIS (West Antarctic Ice Sheet) Divide ice core, Jones et al. (2017) invoked excess diffusion as one of several explanations why diffusion lengths derived from δD signals at  15–18 ka BP exceeded modelled diffusion lengths based on Ramseier's diffusivity by up to 1.6 times. For the high-resolution δD record of the EPICA (European Project for Ice Coring in Antarctica) Dome C core, Pol et al. (2010) regarded vein-driven enhancement of isotopic diffusion to be the cause of strong smoothing and near absence of sub-millennial scale signals across Marine Isotope Stage 19 (MIS 19) – the oldest interglacial identified in that core, at  780 ka BP – where spectrally derived diffusion lengths ( 40 to 60 cm) are several times higher than predicted by the monocrystal diffusivity. Excess diffusion is expected to impact the preservation of deep isotopic signals in the ice cores to be retrieved at Little Dome C, Antarctica, by the European Beyond EPICA project and the Australian Million Year Ice Core project, which aim to obtain records reaching back  1–1.5 Ma and covering the Mid-Pleistocene Transition.

Herein, we revisit Nye's (1998) and Rempel and Wettlaufer's (2003) formulation for vein-mediated excess diffusion, asking “what if the vein water isn't stagnant, but percolates?” In our modelling in Sects. 2 and 3, we show that vein-water flow distorts the isotopic concentrations in ice grains and modifies their isotope exchange with the veins, so that excess diffusion acting on isotopic signals is always increased from the no-flow case. The mechanism causes signals to move relative to the ice, although the age offset of displaced signals is much less than their absolute age. Depending on the water flow velocity, the decay-rate enhancement for decimetre-scale signals can vary by a factor of several to a few hundred – between Rempel and Wettlaufer's and Nye's predictions. This modulation highlights the vein hydrology of ice sheets as a major knowledge gap. In Sect. 4 we explore its impact on signal smoothing at ice-core sites by embedding it in diffusion-length simulations for the GRIP and EPICA ice cores. We show that excess diffusion can undermine diffusion-based temperature reconstructions and the spectral derivation of diffusion lengths from isotope records. We conclude with broader perspectives in Sect. 5.

2 Mathematical model

As in Nye's (1998) and Rempel and Wettlaufer's (2003) studies, our modelling in this section focusses on the interactions between crystals and veins at the grain scale, ignoring the effect of ice deformation on isotopic signals and ignoring diffusion along grain boundaries. Vertical compression will be accounted for in Sect. 4. Background about the diffusion-length theory will be given there.

Figure 1Model geometry and symbols. (a) Ice annular cylinder surrounding a vein. Colour image shows a radial cross-section of the isotopic deviation δ in the ice, exemplifying the signals studied in Figs. 3 to 5. The pattern is distorted by the boundary condition at the vein wall due to isotope advection by vein-water flow. (b) Depth profile of the mean isotopic signal.


First we extend their equations to incorporate vein-water flow. We adopt their idealised geometrical set-up (Fig. 1), which represents ice crystal grains as a vertical annular cylinder, with outer radius b, and the water vein as a hole at its centre, with the vein wall located at the inner radius, r=a (10-6 m); r is the radial coordinate. We distinguish the vein radius a from the radius of curvature of real (convex) vein walls (Nye, 1989; Mader, 1992a; Ng, 2021). In plan view, each cylinder is meant to approximate a unit cell of polycrystalline ice around a vein, so b is taken as the mean grain radius ( 10−3 m). As the original studies assumed, the vein water is kept liquid by a high concentration of dissolved ionic impurities, which lowers the eutectic temperature; horizontal (near-horizontal) veins are disregarded as they cause no (negligible) short-circuiting of depth-varying isotope signals. With the z coordinate axis pointing down, and t denoting time, the concentrations of a trace isotope (18O or D) in the ice grains and in the vein water – Ns(r, z, t) and Nv(z, t), respectively – satisfy the conservation equations


where Ds and Dv are molecular diffusivities in the solid (single crystal) and water. Following the original studies, Nv is assumed independent of r, and we specify Ns/r=0 at r=b as a boundary condition. In Eq. (2), which couples Ns and Nv, the Ds term represents isotope transfer between ice and vein. The final term – our addition to the model – describes advection of Nv by vein water flowing at velocity w (mean value across vein, defined positive downward).

Equilibrium fractionation at the ice–water interface implies αNv/Nv0=Ns|r=a/Ns0, where α ( 1) is the fractionation coefficient, and Nv0 and Ns0 (assumed constant) are the number densities of the major isotope (16O or H) in water and ice. Following the procedure of Rempel and Wettlaufer (2003), we rewrite Nv in Eq. (2) in terms of Ns|r=a, assuming Nv0Ns0, and express Ns as the isotopic deviation δ=δ(r,z,t)=Ns/Ns0-1, thus deriving

(3) δ t = D s 1 r r r δ r + 2 δ z 2 ,

with the boundary conditions

(4) δ r r = b = 0


(5) 2 δ r 2 r = a + 1 - 2 α a δ r r = a - β 2 δ z 2 r = a + w D s δ z r = a = 0 .

The second boundary condition has been derived by eliminating the time derivatives between Eqs. (1) and (2). The dimensionless parameter

(6) β = D v D s - 1 ( > 0 )

quantifies the diffusivity contrast of water to ice. The diffusivities Ds and Dv vary strongly with temperature T, and typically β∼106 (Fig. 2). As detailed in Appendix A, we use Ramseier's (1967) formula for Ds(T), and for Dv(T) we use an extension of Gillen et al.'s (1972) formula that is valid down to 60 C.

Equations (3) to (5) form a partial differential equation model for δ(r,z,t). Rempel and Wettlaufer (2003) assumed 1–2α-1 in Eq. (5), so their model is independent of the fractionation coefficient α and applies equally to δ18O and δD. This is a good approximation because α(18O /16O)  1.0029 and α(D / H)  1.021 (Lehman and Siegenthaler, 1991; O'Neill, 1968; Árnason, 1969). We make the same approximation in most of Sects. 3 and 4, but not in the present derivation, as we need a general model that observes the precise value of α for an analysis about dual-isotope thermometry at the end of Sect. 4.

Equations (3) to (5) encapsulate the short-circuiting effect and differ from Rempel and Wettlaufer's (2003) model by the w term only. When studying the system without water flow (w=0), these authors explained that Nye's (1998) model corresponds to the limit β→∞, as it supposes liquid diffusion so fast that δ along the vein is constant (i.e. perfect short-circuiting); on the other hand, the model of Johnsen et al. (2000) effectively assumes instantaneous radial diffusion in the grains, so that longitudinal diffusion along the vein and in the ice governs the smoothing of signals. Hence, the Johnsen et al. model predicts an excess diffusion equal to Rempel and Wettlaufer's prediction for slow-varying (long) signals, but it is not strictly an approximation of our model at a limit of a defined parameter here. We shall not compare its predictions against our results.

When vein-water flow occurs (w≠0), we expect advection of δ at the vein boundary to perturb δ in the ice (Fig. 1), causing the isotopic signals there to move also. To study how signals behave, we seek a separable solution of the form

(7) δ ( r , z , t ) = δ 0 + δ 1 F ( r ) exp ( - D s ζ t + i k z z ) ,

where ζ=ζR+iζI is a decay-rate parameter for sinusoidal signals with the wavenumber kz (or wavelength λ=2π/kz), and δ0 and δ1 are arbitrary constants representing the background level and amplitude of signals. The amplitude of signals decays at the rate DsζR, whereas their “baseline” decay rate in ice without veins (due to solid diffusion alone) would be Dskz2. Following Nye (1998) and Rempel and Wettlaufer (2003), we let

(8) ζ R = k z 2 + k r 2 = f k z 2 ,

in which the “enhancement factor”

(9) f = 1 + k r 2 k z 2

measures how much faster signals decay in the presence of veins, or equivalently, how much the veins increase the effective diffusivity of the system above Ds. Owing to their short-circuiting effect, excess diffusion operates (f>1) even when w=0. Our main interest is how f varies with w. Note that f, kr, ζR, and ζI are functions of kz.

Figure 2Isotopic diffusivities and their temperature dependence. (a) Arrhenius plot of the self-diffusivity of water Dv. Symbols plot values based on laboratory measurements. Blue curve: our composite exponential in Eq. (A1), fitted to the data of Xu et al. (2016) and used to calculate Dv in this paper. Red curve: quadratic fit by Rempel and Wettlaufer (2003) to the data of Gillen et al. (1972); the dashed portion evaluates their quadratic at T below 31 C, outside its region of applicability. (b) Self-diffusivity of monocrystalline ice Ds, calculated with Ramseier's (1967) empirical formula in Eq. (A2). (c) The liquid-to-solid diffusivity contrast β (=Dv/Ds-1).


Hitherto we seem to be mostly retracing the steps of Nye (1998) and Rempel and Wettlaufer (2003). But a key difference herein – and what distinguishes our findings – is that the decay-rate parameter ζ and the amplitude function F in Eq. (7) are complex numbers when w≠0, since the problem is then no longer symmetric in z. Particularly, a non-zero ζI implies signal migration at the velocity v=ζIDs/kz, and we anticipate F(r)=FR(r)+iFI(r), with the signal phase given by θ(r)=tan-1(FI/FR), varying with radius under the advection. Symmetry considerations for how the system behaves when the vein-water flow direction is reversed predict ζR (hence f) and ζI (hence v) to be even and odd functions of w, respectively.

Now, substituting Eq. (7) into Eqs. (3), (4), and (5) leads to the Bessel equation

(10) F ′′ + F r + ( k r 2 + i ζ I ) F = 0 ,

with the boundary conditions


Suppose kr2+iζIs2=(sR+isI)2, such that

(13) k r 2 = s R 2 - s I 2 and ζ I = 2 s R s I .

Then, in terms of Bessel functions, the analytic solution of Eqs. (10) to (12) is

(14) F ( r ) = J 0 ( s r ) - J 1 ( s b ) Y 1 ( s b ) Y 0 ( s r )

(or any constant multiples), in which s at each wavenumber kz satisfies

(15) 2 α s a s 2 - β k z 2 - i k z w / D s - J 0 ( s a ) Y 1 ( s b ) - Y 0 ( s a ) J 1 ( s b ) J 1 ( s a ) Y 1 ( s b ) - Y 1 ( s a ) J 1 ( s b ) = 0 .

These results are equivalent to those of Rempel and Wettlaufer (2003) when α=1 and w=0. When w≠0, s is complex, and the Bessel functions take complex values.2 For each wavelength λ (or kz), we solve for s of the slowest-decaying mode (with the smallest kr; Nye, 1998) numerically by Newton's Method, taking the left side of Eq. (15) as the function whose root is sought. We then compute kr, ζR, ζI, f and F(r) via Eqs. (13), (8), (9), and (14). The numerical code of these solution steps is given at the online repository of the paper.

Ice-core measurements of δ are often made on horizontal layers spanning multiple ice grains, so it is useful to consider the mean value of δ at each depth in the model (Fig. 1):

(16) δ ¯ ( z , t ) = 1 π b 2 a b 2 π r δ d r + π a 2 δ r = a δ 0 + 2 δ 1 b 2 exp [ - D s ζ t + i k z z ] a b r F ( r ) d r .

This expression shows that the section-mean signal at each wavelength is itself sinusoidal, with the same decay rate, decay-rate enhancement factor, and migration velocity as for the component signals at different radii. The approximation (based on ab) recognises a negligible contribution to the mean signal from the vein water.

3 Results and analysis: excess diffusion at the grain scale

We proceed to analyse computed results to understand the impact of vein-water flow on signal evolution. Notably, we show that advection perturbs δ in such a way that it amplifies the short-circuiting to accelerate signal smoothing, raising f above Rempel and Wettlaufer's (2003) enhancement factor. Through successive numerical experiments, we elucidate the mechanism and key controls on f. We explain relevant properties of the model along the way.

Our experiments here explore signal wavelengths λ in the range 0.001–0.15 m and vein-flow velocities w up to  102 m yr−1, for T=-32C or 52 C. These temperatures resemble those measured in the upper ice column at ice-core sites in central Greenland and central East Antarctica, respectively (e.g. Fig. 8). The higher temperature is close to 241 K, which Rempel and Wettlaufer (2003) chose based on the GRIP site conditions for their calculations. The qualitative dependence of f on the vein and grain radii found by these authors (f increases with a and decreases with b) is unchanged in our model and is not the focus of our study, so we assume constant radii a=1µm and b=1 mm in the experiments. We assume α≡1, the approximation used by Rempel and Wettlaufer (2003), so the results are applicable for either δ18O or δD.

Our range for w is informed by the theoretical estimates of Nye and Frank (1973) for glaciostatically driven water flow through a vein network, listed in their Table 2. We bear in mind that their highest vein-flow velocity (900 m yr−1) assumes a high liquid fraction or porosity ( 10−3) that is probably uncommon for polar ice, so we explore up to values of w an order of magnitude smaller. In the Discussion we comment more on the lack of observations of vein-water flow.

It is noteworthy that Rempel and Wettlaufer (2003) ignored isotope advection by vein-water flow on the basis that the Péclet number Pe=w/Dvkz is small for the signal wavelengths of interest ( dm or less for annual layers). Pe measures the ratio of the advection (fourth) term to the diffusion (third) term in Eq. (5). Both small and large Pe feature in our experiments. As will be seen shortly, although the δ field tends to be modified only in a thin zone by the vein when w≠0, this change can increase f significantly.

Figure 3Patterns of δ calculated with Eqs. (7) to (15) for λ=0.02 m, T=-32C, and w equal to (a) 0, (b) 5, and (c) 50 m yr−1. The enhancement factors in these experiments are f=2.11, 3.24, and 4.26, respectively. Each colour map samples a radial cross-section of the three-dimensional ice annulus in Fig. 1 and is shown with a horizontal exaggeration of 50. Dashed boxes expand on the details near r=0. The vein boundary lies at r=a=1µm. The panel under each map plots the real and imaginary parts of the amplitude function F(r). The panel left of each map plots the δ variations at the vein (red) and in the farthest part of the grain interior (r=b, black).


Figure 3 shows the computed pattern of δ in the ice annulus in three experiments with w=0, 5, and 50 m yr−1 at λ=2 cm and T=-32C. All three signals decay with time; those in Fig. 3b and c migrate downward at constant velocity. We focus on examining the spatial part of the solution in Eq. (7) – the colour maps plot Re[F(r)exp(ikzz)] (or |F(r)|cos(kzz+θ(r)), where θ is the signal phase defined earlier) without the time element. When plotting each pattern, we scale its amplitude such that |F(b)|=1. This choice facilitates comparison of the δ variations along the vein with those at r=b.

Nye's short-circuiting occurs when w=0 (Fig. 3a), as expected. As we move from the grain interior towards the vein, the sinusoidal signals decrease in amplitude sharply following F(r), very close to the vein (dashed box, Fig. 3a). Consequently, isotopes diffuse from the grain interior towards the vein along the peak ridges of the signal, and in the opposite direction along troughs. This pattern of transverse isotopic exchange between vein and ice is driven by fast diffusion along the vein smoothing the signal there and is what causes the entire signal to decay faster than the baseline decay rate, Dskz2 (as governed only by isotope diffusion vertically between the ridges and troughs). The enhancement factor in this experiment is f=2.11, as calculated by Rempel and Wettlaufer (2003).

When we switch on vein-water flow (Fig. 3b, c), f is amplified by a novel effect. Advection shifts the vein signal down relative to the interior, inducing a sheared pattern of δ variations in a thin layer next to the vein boundary, of negative phase (FI<0). At w=5 m yr−1 (Fig. 3b), the sinusoidal variations in the layer have a “tail-like” appearance in colour. Because their phase shift increases towards the vein, there are high lateral gradients in δ at the vein end of the signal peaks and troughs, and they drive a stronger diffusive isotopic exchange between vein and ice (than in the no-flow case), which accelerates the signal decay: f=3.24 in this case. The profile of FR(r) and signal amplitude at r=a are correspondingly reduced.

When w is raised to 50 m yr−1 (Fig. 3c), the shear layer becomes “sheet-like”. Strong advection causes δ outside the vein to interact with δ in the vein much higher up, and the coupling of this with diffusion in the ice diminishes the isotopic variations along and immediately outside the vein to near zero. One can think of the pattern in the layer now as due to the tails of high (low) δ extending far down to cover the next trough (ridge), with δ averaging out sideways by diffusion. Compared to the last experiment, the FR profile is drawn down even more, the transverse isotopic exchange still stronger. The enhancement f=4.26 is nearly maximised, since the signal pattern is close to what it would be at the w→∞ limit, with no variations along the vein as in the Nye model (we confirmed this in experiments that took w above 50 m yr−1).

The last finding implies that, at any signal wavelength, the high flow limit (w→∞) yields the same enhancement as the Nye model limit (β→∞). This equivalence arises because in Eq. (5) w→∞ drives δ/z|r=a to zero, whereas β→∞ drives 2δ/z2|r=a to zero, and both yield the constant vein boundary condition δδ0 to precondition the same isotopic pattern. It follows that f at high flow in our model asymptotically reaches Nye's enhancement factor, and since f at w=0 is minimum (and equal to Rempel and Wettlaufer's enhancement factor), f must take an intermediate value in 0<w<.

Next we consider wavelength control. At fixed w, longer signals experience a higher decay-rate enhancement (though remember their baseline decay rate is lower). Figure 4 demonstrates this with a modified set of experiments, for λ=8 cm, at the same temperature and flow velocities as before. The shear-layer mechanism again operates when w>0. The tails lengthen as w is increased, although the pattern is still in the “tail regime” at 50 m yr−1; at this longer wavelength, a higher w is needed to shift the variations far enough for transition to the “sheet regime”. In all three cases, f is higher (respectively  1.3, 3, and 12 times greater) than before. The reason lies in the relative contribution of (i) vertical diffusion between signal ridges and troughs and (ii) lateral diffusive exchange between vein and ice in driving the signal decay. For signals that are short compared to the grain radius (λb), vertical diffusion dominates over lateral exchange, so vein-water flow increases f minimally via the shear-layer mechanism. For long signals (λb), the lateral exchange is more significant, so the shear-layer mechanism amplifies f more strongly. Note that all three FR profiles in Fig. 4 curve down less than those in Fig. 3, but the attendant reduced lateral exchange rates are still higher than the vertical diffusion rates, which are 16 times less at λ=8 cm than at λ=2 cm. Rempel and Wettlaufer (2003) made similar arguments to explain the wavelength control on f in the system without vein flow. Here we have added the shear-layer mechanism to the considerations.

Figure 4Patterns of δ computed for λ=0.08 m, T=-32C, and w equal to (a) 0, (b) 5, and (c) 50 m yr−1. The enhancement factors in these experiments are f=2.63, 9.01, and 50.2, respectively. The figure has a similar layout as Fig. 3. The inset in panel (a) expands on the variations of FR.


Figure 5Patterns of δ computed for T=-52C at (a) λ=0.02 m, w=0.5 m yr−1, (b) λ=0.02 m, w=5 m yr−1, (c) λ=0.08 m, w=0.5 m yr−1, and (d) λ=0.08 m, w=5 m yr−1. The enhancement factors are f=3.68, 4.27, 14.7, and 51.9, respectively.


Results for colder ice (T=-52C, Fig. 5) predict higher enhancements at the same values of λ and w and shear-layer transitions at lower vein-flow velocities. For both the 2 and 8 cm signals, the tail-to-sheet transition is now largely complete when w reaches 5 m yr−1 (Fig. 5; cf. Figs. 3 and 4). These changes are not due mainly to the change in diffusivity contrast β(=Dv/Ds-1) but rather to the reduced Ds at low temperature (Fig. 2), which raises the importance of vein-flow-assisted lateral isotope exchange compared to vertical diffusion in the grains in smoothing the signal. We study the temperature control more below, where it will be seen that the dominance of these factors Ds and β reverses at low vein-flow velocities.

The mechanism detailed here – initiation of the shear layer by vein-water flow, its progression through the tail and sheet regimes as the magnitude of w is increased, and how the layer isotopic gradients amplify the short-circuiting to accelerate signal decay – is universal across our experiments. To help readers visualise the evolution, we show in Movies S1 and S2 continuous versions of Figs. 3 and 4, with w changing in small steps, covering upward as well as downward water flow.

Having explored the spatial interactions behind the decay-rate enhancement amplification, we report the influences of wavelength and vein-flow velocity more comprehensively by computing curves of f(λ) at fixed w (Fig. 6) and surfaces of f over the λw parameter space (Fig. 7). We do this for T=-32C and T=-52C, plotting log10f and the signal migration velocity v also. Figures 6 and 7 confirm that f increases monotonically with λ and |w|, and, at each λ, f increases from Rempel and Wettlaufer's f value at w=0 towards a maximum (Nye's enhancement factor) as |w|. Importantly, while for centimetre- to decimetre-scale signals f is a few times without vein-water flow, it increases to  101–102 with vein-water flow. The increase is steepest at |w|10 to 20 m yr−1 at 32 C and |w|1 m yr−1 at −52C. Accordingly, for the upper parts of ice cores from central Greenland, West Antarctica and coastal Antarctica, where T-20 to 30 C is common, the extra enhancement above Rempel and Wettlaufer's f is limited until w exceeds a few metres per year (Fig. 6a–b). For the upper parts of ice cores in central East Antarctica, such as at the EPICA Dome C, Dome Fuji, and Vostok ice-core sites, where T-50C, the extra enhancement is already significant at w∼0.5 m yr−1 (Fig. 6d–e). Results computed using the fractionation coefficients for oxygen and deuterium instead of α=1 (Figs. S1 and S2 in the Supplement) differ minimally from those in Figs. 6 and 7.

Figure 6Computed curves of signal decay-rate enhancement factor f, log10f, and signal migration velocity v versus signal wavelength λ, at (a–c) T=-32C and (d–f) T=-52C, for different vein-flow velocities w (labels by the curves, in m yr−1). The lowest curves in f and log10f portray Rempel and Wettlaufer's (2003) enhancement factor. The highest curves approach Nye's (1998) enhancement factor.


The surfaces of f at the two temperatures have similar shapes but different scales in w (Fig. 7). The surface for 52 C is in fact almost exactly a compressed version in the w direction of the surface for 32 C. Movie S3 illustrates this “compressional scaling”, as T varies from 20 to 60 C. The surface always approaches the same profile f(λ) as |w| (also see Fig. 6a, d) to yield Nye's enhancement factor, which does not depend on temperature, because Dv and Ds do not enter the model to influence ζ when w→∞ or β →∞ (Eqs. 3 to 5) (remember the baseline decay rate does increase with T via Ds). However, the surface evolves not merely by compressional scaling. A subtle change also occurs in the valley near w≈0: there, f at 32 C exceeds f at 52 C slightly (the lowest curves in Fig. 6b and e). Consequently, as T is reduced, the surface contracts towards w=0, causing f to rise at moderate to large values of |w| (yielding the earlier result that f decreases with temperature) but to drop near w≈0 (this is too small to be visible in Movie S3).

Figure 7Computed signal decay-rate enhancement factor f, log10f, and signal migration velocity v over the λw parameter space at (a–c) 32 C and (d–f) 52 C. Dashed white curves delineate f(λ) at w=0, which is the enhancement factor in Rempel and Wettlaufer's (2003) model. Circles locate the experiments of Figs. 3, 4, and 5. The curves in Fig. 6 are transects of these surfaces at fixed w.


These temperature controls can be explained by a model scaling analysis. As detailed in Appendix B, with constant vein and grain radii (a and b fixed), three dimensionless parameters govern the signal pattern in the ice annulus and the associated signal decay rate: (i) the ratio of the signal wavelength to the grain radius λ/b, (ii) the ratio of isotope advection by vein-water flow to isotope diffusion in the ice wb/Ds, and (iii) the diffusivity contrast β (the Péclet number considered by Rempel and Wettlaufer, 2003, is a combination of these parameters). It follows that the enhancement factor f has the functional form f(λ/b,wb/Ds,β), whose shape is portrayed by the surfaces in Fig. 7. The influences of λ and w in the first two arguments of this function were explored in earlier experiments. A temperature change affects f via both its second and third arguments, because Ds and β vary with T (Fig. 2). The compressional scaling stems from the second argument, wb/Ds. In contrast, the influence on f by the third argument β (over its range of interest,  106) is weak and prevails only when f is small – near w≈0. There, f decreases when T is reduced from 32 to 52 C because a decrease in β (Fig. 2c) weakens the short-circuiting.

Turning to the migration velocity v (Figs. 6c, f; 7c, f), the model predicts signals to move in the direction of vein-water flow at speeds that reach a maximum at intermediate w and are higher for long signals and at high temperatures, of up to  1 cm kyr−1. Most speeds on the parameter space are much lower. Hence, signal migration is slow, in the sense that long (at least millennial) timescales are needed to displace centimetre- and decimetre-scale annual signals against the ice and other ice-core proxies by a wavelength or more. The relative inaccuracy caused by this on the age scales determined through counting of δ cycles on isotope records is negligible. Compressional scaling applies also to v (Fig. 7, Movie S3), which has the form (Ds/b)g(λ/b,wb/Ds,β) where the function g differs from f (Appendix B). The prefactor Ds/b explains why migration slows as temperature is reduced.

In their calculations, Rempel and Wettlaufer (2003) and Johnsen et al. (2000) accounted for the misorientation of veins from the vertical in the three-dimensional vein network by reducing Dv by a bulk tortuosity factor τf=3. Doing this in our experiments would lower β by ≈3 times and alter the results numerically but not change our qualitative findings.

4 Implications for diffusion-length studies

To explore how much excess diffusion modulated by vein-water flow impacts signal smoothing down the ice column, we simulate diffusion-length profiles for ice-core sites in Greenland and Antarctica, in a forward model testing w. We compare the results against profiles modelled without excess diffusion and query whether they match the level of excess diffusion inferred from isotope records. We also consider the impact on diffusion-length-based temperature reconstructions.

We use the well-established theory of Johnsen (1977) for these calculations, treating what happens below the firn only. In a moving coordinate system where z measures depth below a material horizon in the ice as it descends towards the bed, isotopic signals evolve according to

(17) δ ¯ t = D ( t ) 2 δ ¯ z 2 - ε ˙ z ( t ) z δ ¯ z ,

where t is the age of the horizon, ε˙z (< 0) is the local vertical strain rate, and δ¯ is the section-mean signal in Eq. (16). Since the enhancement factor f applies to this signal (Sect. 2), the bulk-ice isotopic diffusivity in Eq. (17) is given by

(18) D ( t ) = D s ( T ) f ( λ , w ) ,

where the dependence on age arises through T, λ, and other controls of f that vary as the signal descends (β, potentially also w, a, and b). We model a steady-state ice column with fixed temperature and ice velocity profiles. Given the ice thickness, the surface accumulation rate, and the strain-rate profile, the age–depth scale is determined, and signals with the wavelength λ0 at the firn transition shorten to the wavelength λ=λ0S(t)/S(t0) at age t, where S(t)=exp0tε˙z(τ)dτ 1 is the thinning function, and t0 is the firn transition age. The normalisation of S by S(t0), absent in studies that track signals from the ice-sheet surface, accounts for the minor thinning that has taken place by t=t0.

According to Johnsen's (1977) solution of Eq. (17), one can track the amplitudes of different harmonics (Fourier components) of the signal separately3 by using the diffusion length σ, which measures the root mean square distance travelled by diffusing isotopes. Specifically, the squared diffusion length σ2 obeys the differential equation

(19) d σ 2 d t - 2 ε ˙ z ( t ) σ 2 = 2 D ( t ) ,

and each harmonic attenuates by the ratio R=exp(-2π2σ2/λ2) as σ and λ evolve down column. It follows that a harmonic is attenuated strongly when its wavelength shortens to less than σ (e.g. Gkinis et al., 2021). In our simulations, we specify an initial value σ=σfirn at the firn transition – taken from studies of firn isotope diffusion – to circumvent the need to model firn processes. The ratio tracking signal amplitude in the ice is then

(20) R i = exp - 2 π 2 ( σ 2 / λ 2 - σ firn 2 / λ 0 2 ) .

Following Gkinis et al. (2014), we decompose σ2 into a part due to isotopic diffusion in ice, σice2, and another part inherited from firn isotopic diffusion that thins under the compression; thus,

(21) σ 2 ( t ) = σ ice 2 ( t ) + σ firn 2 S ( t ) S ( t 0 ) 2 .

Substituting for σ2 in Eq. (19), using Eq. (18) for D, yields

(22) d σ ice 2 d t - 2 ε ˙ z σ ice 2 = 2 D s f ( λ , w ) ,

with σice=0 at t=t0. This equation is straightforward to solve analytically by an integrating factor (Gkinis et al., 2014) or numerically by the finite-difference method (as done in our simulations below).

In reconstructions of firn temperature (e.g. Gkinis et al., 2014; Holme et al., 2018), Eq. (21) is rearranged as

(23) σ firn 2 = σ obs 2 - σ ice 2 [ S / S ( t 0 ) ] 2 ,

to allow σfirn at a given age to be found from (i) the thinning S, (ii) the modelled value of σice (from Eq. 22), and (iii) the diffusion length σobs estimated from isotopic signals in the ice core – at the same age. Firn isotope diffusion modelling is then used to invert σfirn for temperature. The estimation of σobs involves fitting P0(k)=P0R2=P0exp(-k2σ2) (k=1/λ is the wavenumber) to the power spectral density (PSD) graph of the measured signals, assuming a white-noise input signal at the surface, i.e. constant P0; see Kahle et al. (2018) for different approaches to this estimation.

Equation (22) reveals a notable consequence of excess diffusion (f>1) for the diffusion-length estimation. Besides raising σice (thus σ) to accelerate signal decay, excess diffusion makes σice wavelength dependent: this does not arise if the bulk ice has Ramseier's diffusivity (f≡1), as assumed in most past studies. With excess diffusion, σice varies with λ and the initial wavelength λ0, so the harmonic components have different diffusion-length histories. Their spectral power now decays as exp(k2σ2), where σ decreases with k (rather than being constant), as f increases with λ (Sect. 3). Since f=1 at zero λ only (Figs. 6 and 7), σ exceeds the σ value for monocrystalline ice at all k<∞. Fitting of the resulting non-parabolic PSD thus overestimates σobs, in the context of firn-temperature reconstructions assuming Ramseier's diffusivity for the ice in Eqs. (22) and (23).

More precisely, our model implies that fitting exp(k2σ2) to the PSD decay of the signal to find σ is no longer appropriate when excess diffusion operates; strictly speaking, the fit should be made with exp[-k2(σice2+(σfirnS/S(t0))2)] to find σfirn, with σice (solution of Eq. 22) varying with k and w. We demonstrate the wavelength dependence of σice in simulations below. The dependence does not arise in the firn, because σfirn is not a function of λ, according to models of firn isotope diffusion (Whillans and Grootes, 1985; Johnsen et al., 2000; Gkinis et al., 2014, 2021). Note that when the diffusion length σ develops wavelength dependence, it refers to the root mean square displacement of isotopes only if the signal (at a given depth) has a single wavelength, not if it is composed of multiple harmonics.

We proceed to examine diffusion-length profiles computed for ice-column conditions based on the GRIP site in Greenland and the EPICA Dome C site in Antarctica (Fig. 8). For these sites we specify σfirn=8 and 7 cm, respectively (e.g. Fig. S2 of Gkinis et al., 2014). Most of our model runs assume a=1µm and b=2 mm, although some prescribe a variable grain-radius profile from measurements (Fig. 8d, h). We set α≡1 again, and we set 65 m as the firn transition depth. The age–depth scales yield t0=286.3 a for GRIP and 2872.9 a for EPICA. At the firn transition, the wavelength of annual signals is given approximately by the ice-equivalent surface accumulation rate: 0.23 m yr−1 at GRIP and 0.023 m yr−1 at EPICA.

Figure 8Depth profiles of glaciological variables used in our diffusion-length modelling for the (a–d) GRIP and (e–h) EPICA Dome C ice-core sites. (a, e) Age–depth scale. (b, f) Strain rate ε˙z and the associated thinning function S. (c, g) Ice temperature. (d, h) Grain-radius data (circles) and spline curves used in our modelling. Ice flow at the GRIP site assumes the Dansgaard–Johnsen model with ice thickness H=3029 m, kink height at 1000 m, and surface accumulation rate as=0.23 m yr−1 ice equivalent. Ice flow at the EPICA site assumes the ice submergence velocity wi=mb+(as-mb)(h/H)1.7, where H=3275 m, as=0.023 m yr−1 ice equivalent, h is height above the bed, and mb=0.0008 m yr−1 is the basal melt rate. See Ng (2021) for additional details about these profiles.


Figure 9a–d presents σ profiles simulated at GRIP for the annual signal (λ0=0.23 m) alongside profiles of the enhancement f, ice diffusion length σice, and the ratio Ri tracking signal amplitude. At each depth, σ increases with w via its modulation on σice. Whereas excess diffusion with w=0 (no vein-water flow) raises σ and σice slightly above their values based on Ramseier's diffusivity (f=1), w from several to tens of metres per year increases and modulates σice significantly, with a strong impact on signal decay, down to ≈2300 m depth. From about half way down the column, the amount of excess diffusion diminishes rapidly with depth (f→1) due to severe shortening of the signal and increasing temperature (Sect. 3). The σ and σice profiles converge on the profiles for f=1 near the base of the column because of this, because the firn part of σ is vanishing, and because the long time spent by deep ice at similar strain rate and temperature allows the effects of isotopic diffusion and vertical compression to balance, with σice equilibrating to (-Ds(T)/ε˙z)1/2 in Eq. (22).

Figure 9Computed depth profiles for the GRIP core site of (a, e) the enhancement factor f (representing excess diffusion above the monocrystal diffusivity), (b, f) diffusion length σ, (c, g) ice diffusion length σice, and (d, h) signal amplitude ratio Ri. The number by each curve indicates the vein-water flow velocity w in metres per year (m yr−1). In panels (a) to (d), all model runs study the annual signal (λ0=0.23 m), assuming the grain radius b=2 mm; the dashed curves show results based on Ramseier's (1967) monocrystal diffusivity, i.e. f=1; and the dotted curve in panel (b) shows the thinned firn diffusion length. In panels (e) to (h), blue curves report results based on the variable grain-radius profile in Fig. 8d; red curves report results computed for signals with λ0 a hundred times longer; and grey curves show selected results from panels (a) to (d) for comparison.


The model run for w=20 m yr−1 in Fig. 9a–d mimics the diffusivity enhancement f∼10–30 found for annual signals in the GRIP Holocene ice by Johnsen et al. (1997, 2000), predicting also the reduced excess diffusion in deeper ice (> 1600 m) dating to the Younger Dryas and the last glacial alluded to in those studies. Although blockage of veins by the high dust content in stadial and glacial ice might explain this reduction (Johnsen et al., 1997, 2000), our model predicts a strong dependence of f on the shortening signal wavelength when w>0 (Fig. 6a, b) that provides an explanation. These considerations are unchanged if the run uses variable grain radius (Fig. 9e–h, blue curves for 20 m yr−1). In contrast, the large enhancement f∼10–30 cannot be reproduced with w=0 (Rempel and Wettlaufer's (2003) model) with constant or variable b (Fig. 9e–h), not unless very large vein radii of  20 to 200 µm are assumed.

The wavelength dependence of σ and σice is apparent from two runs that study a signal with an initial wavelength 100 times longer than the annual signal (red curves in Fig. 9e–h; cf. grey curves). The dependence strengthens with w; we analyse its depth variations later. Even with strong excess diffusion (when w=20 m yr−1), the long signal survives much deeper than the annual signal (to ≈2800 m; Fig. 9h) because its baseline decay rate is 104 times lower.

Turning to the EPICA site, our interest is drawn to long signals because annual signals are too short to survive isotopic diffusion in the firn (e.g. for λ=0.023 m, the firn attenuation is R=exp(-2π2σfirn2/λ2)10-80). Figure 10a–c presents model runs for a millennial-scale signal with λ0=23 m. They show a similar modulation of the f and σ profiles by w as seen in the GRIP runs; however, colder ice in the top half of the EPICA column than at GRIP (Fig. 8g, c) means that these profiles are more sensitive to w (Sect. 3). The modulation at EPICA extends to nearer the base of the ice column because the lower strain rate (which shortens signals more slowly and slows the equilibration in Eq. 22) causes σ to approach the f=1 curve slowly. The millennial signal survives into the deepest ice if w≲30 m yr−1 (Fig. 10c).

What conditions at EPICA can produce the long diffusion lengths σ 40 to 60 cm inferred for ice at the depth of MIS 19 ( 3170 m) and thus strong suppression of millennial signals and near-complete absence of sub-millennial signals there? Pol et al. (2010) surmised excess diffusion as necessary. Like the σ profile they simulated, our result for f=1 yields only σ≈16 cm at that depth (Fig. 10b). The runs with excess diffusion at constant b show that w≈80 to 150 m yr−1 generates enough excess diffusion, but the corresponding σ profiles bulge at ≈1000 m to attenuate the millennial signal strongly mid-column (Fig. 10b, c). Other ways of achieving σ 40 to 60 cm in the deepest ice, with small or no bulge in σ that allows a sizeable signal to survive past  2500 m depth but not to  3100 m, are shown in Fig. 10d–g. They assume w profiles ramping up towards the bed – linearly, parabolically, or linearly/nonlinearly near the base – that all require w≳150 m yr−1 in deep ice. If we further consider that Pol et al. (2011) estimated the diffusion length σ 8 cm for ice at MIS 11 ( 395 to 426.7 ka,  2699 to 2799 m) in the EPICA core, then the run using the deep nonlinear w profile (red) best mimics the observations. In any case, excess diffusion unassisted by high vein-water flow at depth cannot explain them. The limited excess diffusion at w= 0, which hardly alters the signal evolution predicted by Ramseier's diffusivity (Fig. 10), is far from sufficient.

Figure 10Computed depth profiles for the EPICA Dome C core site of (a, e) enhancement factor f, (b, f) diffusion length σ, and (c, g) signal amplitude ratio Ri. All model runs investigate a millennial-scale signal with λ0=23 m (1000 × annual layer thickness based on surface accumulation rate). In panels (a) to (c), each run assumes b=2 mm and constant vein-water flow velocity w (value beside each curve, m yr−1). Curves labelled f=1 show results based on Ramseier's diffusivity. Panels (e) to (g) report five model runs able to yield σ of 0.4–0.6 m near the base, assuming the variable w profiles in panel (d). Three runs assume b=2 mm, and two runs assume the variable grain-radius profile in Fig. 8h.


In summary, for both the GRIP and EPICA cores, vein-flow modulation with suitable choices of w can reproduce the levels of excess diffusion inferred from their isotope records. Nye's model (which is approached at high w in Figs. 9 and 10, as w→∞ regains his model; Sect. 3) and Rempel and Wettlaufer's model (w=0) overpredict and underpredict those levels, respectively, in simulations using grain sizes similar to those measured. The reason Rempel and Wettlaufer's model revises down Nye's enhancement factor f was explained in Sects. 2 and 3. Accounting for vein tortuosity (by lowering Dv) weakens the short-circuiting and reduces f further and does not alter these findings.

Figure 11Wavelength dependence of the diffusion length σ (colour, contours in m) at the GRIP site for three parameter settings: (a) w=5 m yr−1, b=2 mm; (b) w=20 m yr−1, b=2 mm; and (c) w=20 m yr−1 with the grain-radius profile in Fig. 8d. In each panel, the maps plot σ versus depth against the initial wavelength of the signal λ0, its initial wavenumber k0, and its wavenumber k at depth. The dashed white lines on the map at far right locate the transects of Fig. 12a.


Figure 12(a) Diffusion length σ as a function of signal wavenumber k, computed at six depths in the GRIP ice column (numbers in colour) in the model runs in Fig. 11c, assuming w= 20 m yr−1 and the variable grain-radius profile. (b) Synthetic power spectral density (PSD) decays based on the curves in panel (a), colour coded as there. Each signal with wavenumber k has k/2π cycles per metre, so the horizontal axes in panels (a) and (b) are equivalent.


And what of the implications for firn temperature reconstructions? Excess diffusion can bias them in two ways, which we discuss here using the GRIP runs as an example.

First, it introduces uncertainty and bias to σobs found by PSD fitting. Figure 11 maps the computed diffusion length versus depth for signals of different initial wavelength λ0 (think of the left-most plot of each panel as showing multiple σ profiles in a stack); in the same way, it maps σ against their initial wavenumber k0=1/λ0 and their wavenumber at depth after accounting for thinning, k=k0S(t0)/S. Whether the simulations use a constant or variable grain radius, σ develops a pronounced wavelength dependence at w=20 m yr−1 (the velocity yielding the level of excess diffusion in GRIP Holocene ice) at most depths, except near the surface and base (Fig. 11b, c). As anticipated, σincreases with λ0 and decreases with k monotonically. Its variations (> 15 % below 700 m; Fig. 12a) mean that the PSD decays [ exp(k2σ2(k))] are not parabolic in k as they might seem (Fig. 12b). Consequently, σobs found by fitting exp(-k2σobs2) to the decays will be misestimated and, as explained earlier, biased too high for existing firn-temperature reconstructions (by an amount dependent on the fitting method).

This issue will arise wherever excess diffusion occurs and in deeper ice below any section with excess diffusion. It will affect reconstructions using the difference between the diffusion lengths of oxygen and deuterium (Simonsen et al., 2011; Holme et al., 2018) as well as those reconstructions using the diffusion length of a single isotope. It can be diagnosed by a statistically significant negative trend in log(PSD)/k2 over k, although real PSD data contain noise and artefacts related to the ice-core isotopic measurements (e.g. Kahle et al., 2018) that may complicate such a test.

Second, excess diffusion affects the calculation of σfirn in Eq. (23) – via the magnitude of σice. Temperature reconstructions based on diffusion length (Gkinis et al., 2014; Holme et al., 2018; Gkinis et al., 2021; Kahle et al., 2021) typically use Ramseier's diffusivity to find σice (i.e. Eq. 22 with f=1), which is justified for ice cores where excess diffusion does not operate. But where it does operate, σice is underestimated, and σfirn and the reconstructed firn temperature are overestimated: accounting for excess diffusion would yield a lower temperature. Put another way, a high σobs at a given depth may result from a large (thinned) σfirn contribution from high firn temperature in the past, from excess diffusion that raised σice, or from both. Robust inversion for σfirn must therefore ascertain the amount of excess diffusion.

The impact on σfirn in single-isotope reconstructions can be gauged using the GRIP runs. Figure 9c and g show that the underestimation of σice by Ramseier's diffusivity is minimal at w=0, but it reaches 0.01 to 0.02 m (to  2000 m depth) if w=20 m yr−1. To gauge how much σfirn is perturbed, we calculate σfirn in Eq. (23) by using σ profiles simulated in forward runs at different w as the input for σobs, but using the profile of σice at f=1 as the other input. These tests thus imitate Ramseier-based reconstructions applied to isotope records influenced by excess diffusion, although they are not comprehensive, as our runs used a fixed imposed history of σfirn≡0.08 m rather than time-varying σfirn. The results show that σfirn is overestimated by an amount increasing with depth and w (Fig. 13). When w=20 m yr−1, σfirn is too high by ∼0.02 m between 1000 and 2000 m depths (ca. 5.5–17 ka), and the overestimation increases from ∼0.01 to 0.03 m through this depth range. With the typical sensitivity of the firn-temperature inversion for Greenland (Fig. S2 of Gkinis et al., 2014; Fig. 3 of Gkinis et al., 2021), a deviation in σfirn of ±0.02 m around 0.08 m translates to a temperature change of ± 5 C. Hence, the firn temperature in this scenario will be overestimated by Ramseier-based reconstructions by several degrees, increasing (with age) through the Holocene and the period since the Last Glacial Maximum.

Figure 13Firn diffusion length σfirn recovered using Eq. (23), with σobs taken from the diffusion-length profiles simulated in forward runs at the GRIP site (for scenarios with excess diffusion at different water-flow velocities w), and with σice calculated by assuming Ramseier's diffusivity for the bulk ice. The forward runs are based on b=2 mm in panel (a) and based on the measured (variable) grain-size profile in panel (b).


In this connection, Gkinis et al. (2014) reconstructed a temperature history from δ18O in the NGRIP (North Greenland Ice Core Project) ice core (retrieved 325 km NNW of GRIP), which they regarded as  3 to 5 C too warm from  8 to 12 ka when compared to other reconstructions (see their Fig. 6). They down-adjusted their temperature results by modifying the thinning function S to reflect lower past accumulation rates. Although the glaciological conditions at NGRIP and GRIP differ, our results show that adjustments of this size are possible with vein-flow modulated excess diffusion. This is not to say that excess diffusion with w∼20 m yr−1 occurred at NGRIP – we are not aware of reports of excess diffusion for that core. However, we recommend that diffusion-based temperature reconstructions assess whether their results could be biased by excess diffusion, besides considering the choices and uncertainties related to firn modelling (Gkinis et al., 2021) and diffusion-length estimation (Kahle et al., 2018). Our test case with w=20 m yr−1 at GRIP is exemplative. Any bias will depend on the magnitude, distribution, and temporal variations of w.

Dual-isotope reconstructions should be less affected by this problem. These reconstructions exploit the different diffusion lengths of δ18O (or δ17O) and δD in the firn, as caused by the different fractionation coefficients (α) for 18O–16O (or 17O–16O) and D–H; specifically, they use the square differential Δσ2=σ2(oxygen) σ2(deuterium) as the proxy for firn temperature (Simonsen et al., 2011; Holme et al., 2018). If the profiles of σice were identical for oxygen and deuterium, as implied by our model of excess diffusion with α≡1 (and by Ramseier-based models), then neither the size of σice nor the bias on σice matter, because σice2 cancels out in the differencing. But the cancellation is imperfect because oxygen and deuterium differ slightly in their fractionation coefficients. To study the effect, we repeated the GRIP runs by using their α values in Eq. (15) to compute the corresponding σice profiles and the ice part of the differential, Δσice2=σice2(oxygen) σice2(deuterium). When w∼10 to 50 m yr−1, Δσice2 reaches  10−5 m2 mid-column (Fig. S3). Since the observed variations in Δσ2 for central Greenland fall in the range 10−4–10−3 m2 (e.g. Figs. 2 and 3 of Simonsen et al., 2011), the ice contribution to the differential can bias dual-isotope reconstructions slightly where excess diffusion operates.

We have not repeated the foregoing analyses for the EPICA ice column because information about its pattern of excess diffusion, which is limited to the σ estimates for MIS 19 and MIS 11 (Pol et al., 2010, 2011), is less complete than at GRIP.

5 Conclusions and outlook

In this paper, we described a mechanism whereby vein-water flow amplifies the short-circuiting conceived by Nye (1998), enhancing the rate of isotopic diffusion in polycrystalline ice above the rate predicted by Rempel and Wettlaufer's (2003) model. Our simulations demonstrate its profound impact on signal smoothing in ice where the vein-water flow velocity w reaches  101–102 m yr−1. We explained why vein-flow modulated excess diffusion biases the spectral estimation of diffusion lengths from isotope records, as well as diffusion-based palaeothermometry at ice-core sites. Our findings contribute insights essential for developing robust interpretation of ice-core isotope records.

Potential modulation of excess diffusion by vein-water flow means that ice-core isotopic signals may have been altered in more complex ways than previously thought. Where modulation occurs, neither Ramseier's (1967) formula nor Rempel and Wettlaufer's (2003) model describe the bulk-ice isotopic diffusivity, and we caution their use in ice-core analysis.

Our GRIP and EPICA simulations (Sect. 4) represent a detailed exploration of the impact of excess diffusion on diffusion-length profiles and probe the conditions behind the levels of excess diffusion inferred for those cores. Their expository nature should be emphasised: we have not fitted the observations precisely, and the vein-flow velocities in our model runs are only trial values in a sensitivity analysis. Several reasons preclude bespoke modelling for accurate validation or prediction at present. First, the model idealises the system geometry. Based on Nye's and Rempel and Wettlaufer's set-up (Fig. 1), it conceptualises the vein network in three dimensions by using a uniform array of cells and approximates the vein cross-section – which in reality consists of three convex faces (Nye, 1989; Mader, 1992a) – as circular. It also assumes that veins are static features, even though grain-boundary migration implies a continual slow motion of veins relative to crystal matrices (Ng, 2021), which will perturb the isotope concentration fields within ice grains. Second, the model ignores isotope diffusion along grain boundaries (Johnsen et al., 2000), which may further enhance excess diffusion in fine-grained ice (Jones et al., 2017). Third, some material properties serving as model inputs are subject to uncertainty. In particular, the formulae we use for the diffusivities Ds and Dv may only be reliable to within a factor of several; they neglect the potential influence of pressure and dissolved ionic impurities (which may change the fractionation coefficients also). The uncertainty involved is probably not dissimilar in magnitude to that associated with the geometrical approximation, and our qualitative findings (predicated on Dv/Ds1) are robust against it; how numerical results vary in our model if Ds and Dv are altered is also straightforward to compute. However, this uncertainty impacts all key modelling studies of excess diffusion to date (Nye, 1998; Johnsen et al., 2000; Rempel and Wettlaufer, 2003; the present study), so there is a clear need for more comprehensive laboratory determination of both diffusivities. Fourth, in simulations made for specific ice-core sites, two factors besides w are challenging to constrain: (i) suppression of the short-circuiting due to blockage of veins by dust particles and bubbles and (ii) spatial variation in vein radius. The vein size in ice at thermodynamic equilibrium should depend locally on the mean grain size, temperature, and amount of dissolved ionic impurities in the veins (Nye, 1991; Mader, 1992b; Dani et al., 2012). Consequently, factors (i) and (ii) suggest possible influences by changing vein-impurity signals (Ng, 2021) and changing distributions of bubbles and solid particles on excess diffusion and the smoothing of isotopic signals. Extending the model for these controls and the detailed geometry of the vein network as grain size and texture evolve are worthwhile avenues for further research.

A striking realisation from our study is how little is known about the vein-scale hydrology of ice sheets. The vein-flow velocity w is needed for predicting excess diffusion with the model or validating model runs made to match diffusion lengths measured from ice cores. However, reliable prediction of the size and pattern of w at ice-core sites is out of reach. Our only handle on w is Nye and Frank's (1973) theory for glaciostatically driven porewater flow. This theory calculates the rate of (Poiseuille) flow through veins under the hydropotential gradient (ρwρi)g, where ρw is water density, ρi is ice density, and g is gravitational acceleration. It yields a large range of plausible w because the ice porosity – a key input that determines the vein size for a given mean grain size – is uncertain for polar ice. When modelling the profile of w above subglacial lakes, Rempel (2005) combined Nye and Frank's (1973) theory with an equation of vein-equilibrium thermodynamics to constrain the porosity through the influence of dissolved impurities. Although this approach can predict specific values of w, it requires knowledge (or assumption) about the amount of ionic impurities partitioned to the vein network, which is not resolved by most ice-core analytical measurements. Besides, neither Rempel's nor Nye and Frank's model has been observationally tested. A separate theory by Nye (1976) treats vein-water flow in detail but analyses only its stability, not flow rates. Thus, currently, knowledge about w in polar ice is limited to a few embryonic theories, and we cannot evaluate whether the range of w values found to modulate excess diffusion sensitively are common at ice-core sites.

Direct measurements of vein-water flow in polar ice are critical for progress, for informing studies of excess diffusion and hydrological modelling. Vein size measurements at low temperature are also desirable, because the observations of Mader (1992b) (vein widths 10–100 µm) were made only within 1 C below the melting point. It may be possible to measure w by innovating on nuclear magnetic resonance (NMR) and Doppler-based techniques. Laboratory studies should couple vein size and flow measurements with experiments where water percolates through ice containing isotopic signals. In such instances, one can use LA-ICP-MS (laser ablation inductively coupled plasma mass spectrometry; Bohleber et al., 2021) to try to discern the occurrence of excess diffusion, by studying the relationship between isotopic signals in the crystals and those in the veins and by looking for the flow-induced “shear layer” in isotopic concentration near triple junctions predicted in Sect. 3. Hitherto, no experiments have been made to demonstrate Nye's short-circuiting effect, even for the case without water flow.

Because the amplification of excess diffusion is due to isotopic gradients caused by vein-flow advection (Sect. 3), our model implies that any non-zero pattern of w will accelerate the smoothing of isotopic signals. The water percolation need not be vertical or unidirectional, nor occur on long length scales (as experimented herein). We think that recrystallisation in deforming ice will cause nonuniform vein-water flow at the grain scale, although no theories yet address this process quantitatively, and we do not know the flow velocities involved. Notably, polygonisation and strain-induced migration recrystallisation that reconfigure grain boundaries must create new vein segments while eliminating others, thus inducing local water flow superposed on any long-range transport (e.g. the downward percolation envisaged by Nye and Frank, 1973); some local water flow might occur even where veins are blocked. If this hypothesis is correct, then the coupling between recrystallisation processes and isotopic diffusion is more complicated than the irregular shape and motion of existing veins, and the rate and mode of ice deformation will affect the level of excess diffusion.

It is important to establish whether the mechanisms of vein-mediated isotopic diffusion and vein-flow-induced enhancement physically operate as theorised, besides yielding the right levels of excess diffusion. The experiments mentioned above using ice samples with LA-ICP-MS provide a direct test. Tentatively, bearing in mind the difficulty of determining vein-size variations and vein-network connectivity, one might be able to test the theory with ice-core records, by comparing their isotopic signals against impurity signals considered immobile or less mobile, to identify signs of isotopic signal migration (Sect. 3). Portions of ice cores that the theory deems to be affected differently by excess diffusion, such as ice with high and low dust content and/or dissolved impurities, as well as transition regions, might be worthwhile targets for study.

Future studies should analyse the pattern of excess diffusion in multiple ice cores systematically to help unravel its diverse controls; this is not least because the origin of excess diffusion at the GRIP and EPICA sites remains elusive. High-resolution isotopic measurements on ice cores that are becoming more common (e.g. Steig et al., 2021) will aid this effort. Our modelling shows that where excess diffusion occurs, one should be able to quantify and remove the biases on spectrally derived diffusion lengths and firn-temperature reconstructions by calculating σice and its wavelength dependence with Eq. (22). Although this is not possible to do until w can be predicted for each ice-core site, diffusion-based studies should scrutinise their isotope records for signs of excess diffusion by studying the decay rate of signals at specific (unthinned) wavelengths (Johnsen et al., 1997, 2000) and testing for significant trends on log(PSD)/k2 data (Sect. 4).

Appendix A: Diffusivities Dv and Ds

Following Rempel and Wettlaufer (2003), we use a relation for Dv(T) based on experimental measurements on supercooled water, even though strictly speaking the vein water is liquid – in thermodynamic equilibrium with ice – due to the presence of dissolved impurities and interfacial curvature (Mulvaney et al., 1988; Nye, 1991; Ng, 2021), not because of supercooling.

Rempel and Wettlaufer (2003) fitted a quadratic to the self-diffusion coefficients measured by Gillen et al. (1972) with an NMR method for supercooled water down to -31C. At lower T, their quadratic is not meant to apply, especially as it predicts Dv to increase again (Fig. 2a, dashed red curve). One cannot extrapolate the trend of Gillen et al.'s (1972) data to those temperatures, as there has been much debate about a possible liquid–liquid phase transition in the so-called “no man's land” of deeply supercooled water near 228 K (-45C), which is thought to be responsible for various observed anomalies in the properties of water. Reviews of this topic have been given by Amann-Winkel et al. (2016), Handle et al. (2017), and Hestand and Skinner (2018).

Recently, Xu et al. (2016) measured the rate of growth of ice crystals into supercooled water by using a pulsed-laser heating technique and used the results with the Wilson–Frenkel model (Wilson, 1900; Frenkel, 1946) to derive new estimates of Dv, down to 125 K. From their growth and diffusivity data, they inferred no thermodynamic transitions or singularities in no man's land. Herein, we fit their Dv values between 12.8 and 60.8 C (blue circles, Fig. 2a), which exhibit a slope break on the Arrhenius plot at -38C, by using the composite exponential

(A1) D v = 1 1 1.085 × 10 - 6 exp - 1870 T + 1 2.942 × 10 7 exp - 9474 T m 2 s - 1

(blue curve, Fig. 2a). Although their experiment was conducted in an ultra-high vacuum, we use Eq. (A1) without pressure correction in our model, because combined atmospheric and glaciostatic pressures should increase Dv by a few percent only (Dv increases by 10 %–40 % as pressure rises from 0.1 to 100 MPa; Prielmeier et al., 1988) and because most of the data points of Xu et al. actually lie above those of Gillen et al. (Fig. 2a).

For the self-diffusivity of ice Ds, we use Ramseier's (1967) formula

(A2) D s = 9.1 × 10 - 4 exp - 7.2 × 10 3 T m 2 s - 1 ,

which is based on experiments on ice monocrystals down to 35.9 C. We are not aware of measurements of Ds at lower temperatures, where diffusion is too slow to be easily determinable on experimental timescales. However, there is no reason to expect the activation energy of volume self-diffusion ( 60 kJ mol−1 according to Eq. A2) to change drastically at low T. Examples of the use of Eq. (A2) for T as low as -50C appear in earlier studies (e.g. Pol et al., 2010; Grisart et al., 2022).

Appendix B: Scaled model

By non-dimensionalising the independent variables with r*=r/b, z*=z/b, and t*=t/(b2/Ds), Eqs. (3) to (5) become

(B1) δ t = 1 r r r δ r + 2 δ z 2 ,

with the boundary conditions


The parameters

(B4) ε = a b ( 1 ) and γ = w b D s

signify the dimensionless vein radius and the dimensionless vein-water flow velocity, respectively. We have assumed α=1. If we rescale the decay-rate parameter, wavelength, and wavenumbers according to ζ*=b2ζ, λ*=λ/b, and [kz*, kr*] = b[kz,kr] (thus ensuring λ*=2π/kz*), then the solution in Eq. (7) becomes

(B5) δ ( r , z , t ) = δ 0 + δ 1 F ( r ) exp ( - ζ t + i k z z ) ,

and the enhancement factor f=1+(kr*/kz*)2 retains the form in Eq. (9).

This scaled model implies f=f(ε, λ*, γ, β) or f(λ*, γ, β) at fixed ε. The signal migration velocity v (=ζIDs/kz dimensionally) is given by (Ds/b)Im(ζ*)/kz* or (Ds/b)g, where g is another function of the same parameters. Indeed, one can solve the scaled model by the method of Sect. 2 (replacing a by ε, b and Ds by 1, and w by γ) and, after computing f and g, infer the dimensional controls on f and v (considered in Sect. 3 and Fig. 7) through the parameter and scaling definitions.

Code and data availability

The MATLAB code for evaluating the model and the computed data are archived at (Ng, 2023a).

Video supplement

Movies S1–S3 are available at (Ng, 2023b).


Figures S1–S3 and Movies S1–S3 are available at (Ng, 2023b). The supplement related to this article is available online at:

Competing interests

The author has declared that there are no competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Special issue statement

This article is part of the special issue “Ice core science at the three poles (CP/TC inter-journal SI)”. It is not associated with a conference.


I thank Greg Kimmel for providing the diffusivity data from the study by Xu et al. (2016), Andrew Sole and Adam Hepburn for proofreading the submitted paper, and the University of Sheffield for funding contribution towards publication. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Review statement

This paper was edited by Mathieu Casado and reviewed by Kurt Cuffey and one anonymous referee.


Abramowitz, M. and Stegun, I. A.: Handbook of Mathematical Functions with Formulas, Graphs and Mathematical Tables, 10th Edn., National Bureau of Standards, Washington, DC, ISBN 0471800074, 1972. 

Amann-Winkel, K., Böhmer, R., Fujara, F., Gainaru, C., Geil, B., and Loerting, T.: Colloquium: Water's controversial glass transitions, Rev. Mod. Phys., 88, 011002,, 2016. 

Árnason, B.: Equilibrium constant for the fractionation of deuterium between ice and water, J. Phys. Chem., 73, 3491–3494, 1969. 

Augustin, L., Barbante, C., Barnes, P. R. F., et al.: Grain radius from selected samples of the EPICA Dome C ice core EDC, 110-3100 metres, PANGAEA [data set],, 2004. 

Bohleber, P., Roman, M., Šala, M., Delmonte, B., Stenni, B., and Barbante, C.: Two-dimensional impurity imaging in deep Antarctic ice cores: snapshots of three climatic periods and implications for high-resolution signal interpretation, The Cryosphere, 15, 3523–3538,, 2021. 

Casado, M., Münch, T., and Laepple, T.: Climatic information archived in ice cores: impact of intermittency and diffusion on the recorded isotopic signal in Antarctica, Clim. Past, 16, 1581–1598,, 2020. 

Cuffey, K. M. and Steig, E. J.: Isotopic diffusion in polar firn: Implications for interpretation of seasonal climate parameters in ice-core records, with emphasis on central Greenland, J. Glaciol., 44, 273–284,, 1998. 

Dani, K. G., Mader, H. M., Wolff, E. W., and Wadham, J. L.: Modelling the liquid-water vein system within polar ice sheets as a potential microbial habitat, Earth Planet. Sc. Lett., 333–334, 238–249,, 2012. 

Frenkel, J.: Kinetic Theory of Liquids, The Clarendon Press, Oxford, 1946. 

Gillen, K. T., Douglass, D. C., and Hoch, M. J. R.: Self-diffusion in liquid water to 31 C, J. Chem. Phys., 57, 5117–5119, 1972. 

Gkinis, V., Simonsen, S. B., Buchardt, S. L., White, J. W. C., and Vinther, B. M.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years – Glaciological and paleoclimatic implications, Earth Planet. Sc. Lett., 405, 132–141,, 2014. 

Gkinis, V., Holme, C., Kahle, E. C., Stevens, M. C., Steig, E. J., and Vinther, B. M.: Numerical experiments on firn isotope diffusion with the Community Firn Model, J. Glaciol., 67, 450–472,, 2021. 

Grisart, A., Casado, M., Gkinis, V., Vinther, B., Naveau, P., Vrac, M., Laepple, T., Minster, B., Prié, F., Stenni, B., Fourré, E., Steen-Larsen, H. C., Jouzel, J., Werner, M., Pol, K., Masson-Delmotte, V., Hoerhold, M., Popp, T., and Landais, A.: Sub-millennial climate variability from high-resolution water isotopes in the EPICA Dome C ice core, Clim. Past, 18, 2289–2301,, 2022. 

Hammer, C. U., Clausen, H. B., Dansgaard, W., Gundestrup, N., Johnsen, S. J., and Reeh, N.: Dating of Greenland ice cores by flow models, isotopes, volcanic debris, and continental dust, J. Glaciol., 20, 3–26, 1978. 

Handle, P. H., Loerting, T., and Sciortino, F.: Supercooled and glassy water: Metastable liquid(s), amorphous solid(s), and a no-man's land, P. Natl. Acad. Sci. USA, 114, 13336–13344,, 2017. 

Hestand, N. J. and Skinner J. L.: Perspective: Crossing the Widom line in no man's land: Experiments, simulations, and the location of the liquid-liquid critical point in supercooled water, J. Chem. Phys., 149, 140901,, 2018. 

Holme, C., Gkinis, V., and Vinther, B. M.: Molecular diffusion of stable water isotopes in polar firn as a proxy for past temperatures, Geochim. Cosmochim. Ac., 225, 128–145, 2018. 

Johnsen, S. J.: Stable isotope homogenization of polar firn and ice, International Association of Hydrological Sciences Publication 118, Symposium at Grenoble 1975: Isotopes and Impurities in Snow and Ice, 210–219, 1977. 

Johnsen, S. J., Clausen, H. B., Dansgaard, W., Gundestrup, N. S., Hammer, C. U., Andersen, U., Andersen, K. K., Hvidberg, C. S., Dahl-Jensen, D., Steffensen, J. P., Shoji, H., Sveinbjörnsdóttir, Á. E., White, J., Jouzel, J., and Fisher, D.: The δ18O record along the Greenland Ice Core Project deep ice core and the problem of possible Eemian climatic instability, J. Geophys. Res.-Oceans, 102, 26397–26410,, 1997. 

Johnsen, S. J., Clausen, H. B., Cuffey, K. M., Hoffmann, G., Schwander, J., and Creyts, T.: Diffusion of stable isotopes in polar firn and ice: the isotope effect in firn diffusion, in: Hondoh, T., Physics of ice core records. Sapporo, Hokkaido University Press, 121–140, ISBN 4832902822, 2000. 

Jones, T. R., Cuffey, K. M., White, J. W. C., Steig, E. J., Buizert, C., Markle, B. R., McConnell, J. R., and Sigl, M.: Water isotope diffusion in the WAIS Divide ice core during the Holocene and last glacial, J. Geophys. Res.-Earth, 122, 290–309,, 2017. 

Kahle, E. C., Holme, C., Jones, T. R., Gkinis, V., and Steig, E. J.: A generalized approach to estimating diffusion length of stable water isotopes from ice-core data, J. Geophys. Res.-Earth, 123, 2377–2391,, 2018. 

Kahle, E. C., Steig, E. J., Jones, T. R., Fudge, T. J., Koutnik, M. R., Morris, V. A., Vaughn, B. H., Schauer, A. J., Max Stevens, C., Conway, H., Waddington, E. D., Buizert, C., Epifanio, J., and White, J. W. C.: Reconstruction of temperature, accumulation rate, and layer thinning from an ice core at South Pole, using a statistical inverse method, J. Geophys. Res.-Atmos., 126, e2020JD033300,, 2021. 

Küttel, M., Steig, E. J., Ding, Q., Monaghan, A. J., and Battisti, D. S.: Seasonal climate information preserved in West Antarctic ice core water isotopes: relationships to temperature, large-scale circulation, and sea ice, Clim. Dynam., 39, 1841–1857,, 2012. 

Laepple, T., Münch, T., Casado, M., Hoerhold, M., Landais, A., and Kipfstuhl, S.: On the similarity and apparent cycles of isotopic variations in East Antarctic snow pits, The Cryosphere, 12, 169–187,, 2018. 

Lehman, M. and Siegenthaler, U.: Equilibrium oxygen- and hydrogen-isotope fractionation between ice and water, J. Glaciol., 37, 23–26, 1991. 

Mader, H. M.: Observations of the water-vein system in polycrystalline ice, J. Glaciol., 38, 333–347, 1992a. 

Mader, H. M.: The thermal behaviour of the water-vein system in polycrystalline ice, J. Glaciol., 38, 359–374, 1992b. 

Mulvaney, R., Wolff, E. W., and Oates, K.: Sulphuric acid at grain boundaries in Antarctic ice, Nature, 331, 247–249, 1988. 

Ng, F. S. L.: Pervasive diffusion of climate signals recorded in ice-vein ionic impurities, The Cryosphere, 15, 1787–1810,, 2021. 

Ng, F.: Numerical code of the paper “Isotopic diffusion in ice enhanced by vein-water flow”, Figshare [data set and code],, 2023a. 

Ng, F.: Supplement of the paper “Isotopic diffusion in ice enhanced by vein-water flow”, The University of Sheffield [video],, 2023b. 

Nye, J. F.: Water flow in glaciers: jökulhlaups, tunnels and veins, J. Glaciol., 17, 181–207, 1976. 

Nye, J. F.: The geometry of water veins and nodes in polycrystalline ice, J. Glaciol., 35, 17–22, 1989. 

Nye, J. F.: Thermal behaviour of glacier and laboratory ice, J. Glaciol., 37, 401–13, 1991. 

Nye, J. F.: Diffusion of isotopes in the annual layers of ice sheets, J. Glaciol., 44, 467–468, 1998. 

Nye, J. F. and Frank, F. C.: Hydrology of the intergranular veins in a temperate glacier, International Association of Scientific Hydrology Publication 95, Symposium at Cambridge 1969 – Hydrology of Glaciers, 157–161, 1973. 

O'Neill, J. R.: Hydrogen and oxygen isotope fractionation between ice and water, J. Phys. Chem., 72, 3683–3684, 1968. 

Pol, K., Masson-Delmotte, V., Johnsen, S., Bigler, M., Cattani, O., Durand, G., Falourd, S., Jouzel, J., Minster, B., Parrenin, F., Ritz, C., Steen-Larsen, H. C., and Stenni, B.: New MIS 19 EPICA Dome C high resolution deuterium data: Hints for a problematic preservation of climate variability at sub-millennial scale in the “oldest ice”, Earth Planet. Sc. Lett., 298, 95–103,, 2010. 

Pol, K., Debret, M., Masson-Delmotte, V., Capron, E., Cattani, O., Dreyfus, G., Falourd, S., Johnsen, S., Jouzel, J., Landais, A., Minster, B., and Stenni, B.: Links between MIS 11 millennial to sub-millennial climate variability and long term trends as revealed by new high resolution EPICA Dome C deuterium data – A comparison with the Holocene, Clim. Past, 7, 437–450,, 2011. 

Prielmeier, F. X., Lang, E. W., Speedy, R. J., and Lüdemann, H.-D.: The pressure dependence of self diffusion in supercooled light and heavy water, Ber. Bunsenges. Phys. Chem., 92, 1111–1117, 1988.  

Ramseier, R. O.: Self-diffusion of tritium in natural and synthetic ice monocrystals, J. Appl. Phys., 38, 2553–2556, 1967. 

Rempel, A.: Englacial phase changes and intergranular flow above subglacial lakes, Ann. Glaciol., 40, 191–194, 2005. 

Rempel, A. W. and Wettlaufer, J. S.: Isotopic diffusion in polycrystalline ice, J. Glaciol., 49, 397–406, 2003. 

Simonsen, S. B., Johnsen, S. J., Popp, T. J., Vinther, B. M., Gkinis, V., and Steen-Larsen, H. C.: Past surface temperatures at the NorthGRIP drill site from the difference in firn diffusion of water isotopes, Clim. Past, 7, 1327–1335,, 2011. 

Steig, E. J., Jones, T. R., Schauer, A. J., Kahle, E. C., Morris, V. A., Vaughn, B. H., Davidge, L., and White, J. W. C.: Continuous-flow analysis of δ17O, δ18O, and δD of H2O on an ice core from the South Pole, Front. Earth Sci., 9, 640292,, 2021. 

Thorsteinsson, T., Kipfstuhl, J., and Miller, H.: Textures and fabrics in the GRIP ice core, J. Geophys. Res., 102, 26583–26599, 1997. 

Vinther, B. M., Clausen, H. B., Johnsen, S. J., Rasmussen, S. O., Andersen, K. K., Buchardt, S. L., Dahl-Jensen, D., Seierstad, I. K., Siggaard-Andersen, M.-L., Steffensen, J. P., Svensson, A., Olsen, G., and Heinemeier, J.: A synchronized dating of three Greenland ice cores throughout the Holocene, J. Geophys. Res.-Atmos., 111, D13102,, 2006. 

Whillans, I. M. and Grootes, P. M: Isotopic diffusion in cold snow and firn, J. Geophys. Res., 90, 3910–3918,, 1985. 

Wilson, H. A.: On the velocity of solidification and viscosity of supercooled liquids, Lond. Edinb. Dublin Philos. Mag. J. Sci., 50, 238–250, 1900. 

Xu, Y., Petrika, N. G., Smith, R. S., Kay, B. D., and Kimmel, G. A.: Growth rate of crystalline ice and the diffusivity of supercooled water from 126 to 262 K, P. Natl. Acad. Sci. USA, 113, 14921–14925,, 2016. 

Zheng, M., Sjolte, J., Adolphi, F., Vinther, B. M., Steen-Larsen, H. C., Popp, T. J., and Muscheler, R.: Climate information preserved in seasonal water isotope at NEEM: relationships with temperature, circulation and sea ice, Clim. Past, 14, 1067–1078,, 2018. 


Recounting the same analysis, Johnsen et al. (2000) later reported 30 times.


One hopes to express the solution of Eq. (10) in terms of real transcendental functions, as in the case kr>0, ζI=0 (which gives Bessel functions with real arguments) or the case kr=0, ζI<0 (which gives Kelvin functions); see Abramowitz and Stegun (1972). For kr and ζI both non-zero, however, we have not found such functions in the literature and need to evaluate the Bessel functions in Eqs. (14) and (15) for complex arguments. This is done straightforwardly in MATLAB.


This is because their wavelengths follow different histories of shortening. To see this, notice Eq. (17) has the characteristic velocity dz/dt=ε˙zz, which motivates a change of the depth variable to Z=z/S(t). Changing the time variable also via τ=0tD(t)/S(t)2dt converts Eq. (17) to the classic diffusion equation δ/τ=2δ/Z2 , whose Fourier components evolve independently.

Short summary
The stable isotopes of oxygen and hydrogen in ice cores are routinely analysed for the climate signals which they carry. It has long been known that the system of water veins in ice facilitates isotopic diffusion. Here, mathematical modelling shows that water flow in the veins strongly accelerates the diffusion and the decay of climate signals. The process hampers methods using the variations in signal decay with depth to reconstruct past climatic temperature.