Constraints on the delta H-2 diffusion rate in firn from field measurements at Summit, Greenland

. We performed detailed 2 H isotope diffusion measurements in the upper 3 m of ﬁrn at Summit, Greenland. Using a small snow gun, a thin snow layer was formed from 2 H-enriched water over a 6 × 6 m 2 area. We followed the diffusion process, quantiﬁed as the increase of the δ 2 H diffusion length, over a 4-year period, by retrieving the layer once per year by drilling a ﬁrn core and slicing it into 1 cm layers and measuring the δ 2 H signal of these layers. We compared our experimental ﬁndings to calculations based on the model by Johnsen et al. (2000) and found substantial differences. The diffusion length in our experiments increased much less over the years than in the model. We discuss the possible causes for this discrepancy and conclude that several aspects of the diffusion process in ﬁrn are still poorly constrained, in particular the tortuosity.


Introduction
The relative abundance of the stable isotopes 2 H and 18 O in ice cores is one of the most powerful proxies of the palaeotemperature over the last 800 kyr (Jouzel and EPICA community, 2007).The global meteoric water cycle acts as a global-scale isotope distillation system, through a continuous process of evaporation and condensation.It leads to a depletion of the abundance of heavier isotopes in the water molecules, which depletion increases with higher latitudes or rather, in fact, with lower temperature.In the polar regions, especially on Antarctica and Greenland, the precipitation containing this temperature-dependent isotope content is conserved, and by drilling deep ice cores on strategic places on these ice caps the precipitation of over 100 000 years (Dahl-Jensen et al., 2013;Johnsen et al., 2001) or even close to a million years (Oerter et al., 2004;Stenni et al., 2010) can be recovered.Accurate, high-spatial-resolution isotope abundance measurements on this ice core material then reveal the "proxy" temperature signal.
As said, the signals are proxies, implying that their relation with past temperature is solid, but not necessarily linear, and not even of a constant character through time or in space.Most of this proxy character is due to the complex and time-dependent relation between temperature and circulation patterns in the atmosphere, influencing the behaviour of the isotope distillation system.Some of it, however, is due to processes that influence the isotopic abundance pattern after deposition.Apart from processes acting only on very fresh snow still at or close to the surface (firn ventilation and sublimation; Neumann and Waddington, 2004;Sokratov and Golubev, 2009), the main process is diffusion.This smears out, and can eventually wash away, spatial gradients in the isotope abundances.Diffusion takes place in the solid-ice phase, but especially, at a rate some 5 orders of magnitudes higher, in the firn phase.In the latter case, the diffusion process is governed by the continuous evaporation and condensation of ice particles into and from the air channels.While in the vapour phase, the water molecules travel a certain distance before freezing back to the solid phase.The process is relatively fast, especially in the first years after deposition, when the firn density is still low, and summer temperatures still affect the firn temperature.On the central Greenland ice sheet, the Published by Copernicus Publications on behalf of the European Geosciences Union.L. G. van der Wel et al.: Constraints on the δ 2 H diffusion rate in firn from field measurements at Summit isotope diffusion process decreases the seasonal cycle amplitude by typically a factor of 5 and effectively influences all timescales of decadal variability and shorter (Andersen et al., 2006a;Vinther et al., 2006).
Isotope diffusion in firn was discovered by Langway (1967) and is well understood in the qualitative sense.Quantification, however, is more difficult, as the rate depends critically on the density of the firn and, moreover, on the tortuosity of the firn, i.e. the shape and size of the air channels between the ice crystals.Furthermore even grain size might play a role.Several laboratory experiments have been performed on the isotope diffusion rate (Jean-Baptiste et al., 1998;Pohjola et al., 2007;van der Wel et al., 2011a), and expressions for the rates that include a parameterization for the tortuosity dependence have been formulated (Cuffey and Steig, 1998;Johnsen et al., 2000;Whillans and Grootes, 1985) and tested (Pohjola et al., 2007;van der Wel et al., 2011a).
Corrections for isotope diffusion (usually to reconstruct the original precipitation signal: "back diffusion") are currently performed routinely.The last publication that we are aware of that shows the "raw", uncorrected isotope signals together with the corrected ones is by Johnsen (1977).So, in spite of the fact that the process is difficult to quantify, back diffusion is applied routinely for the interpretation of ice cores, even for short(er) timescales (Bolzan and Pohjola, 2000;Vinther et al., 2006).
As the process is difficult to quantify, we found the isotope diffusion rate worth further investigation, especially since (1) isotope diffusion has gained renewed attention, after the discovery (by Johnsen et al., 2000) that the difference of the diffusion rate for 2 H and 18 O ("differential diffusion") is dependent only on the temperature of the firn.Differential diffusion has thus the potential of becoming a powerful new palaeotemperature proxy by itself (Simonsen et al., 2011;van der Wel, 2012).(2) All laboratory experiments on diffusion rates have been performed on artificially produced "firn": shaved ice flakes (Jean-Baptiste et al., 1998;Pohjola et al., 2007) or, at best, snow produced by a snow gun (van der Wel et al., 2011a).However, it was realized that the tortuosity, the 3-D shape of the air channels between the crystals, is important for the diffusion rate, and it is well conceivable that this differs considerably between artificial "snow" and real snow.
(3) Most laboratory experiments have concentrated on the high-density regime, where tortuosity effects are most pronounced but where the diffusion process is slow.Here, we concentrate on the initial phase of firnification where diffusion is fastest.
For these three reasons we decided to design and perform a field experiment in which we could measure the 2 H isotope diffusion rate for real snow.Using a snow gun, we produced a thin layer of artificially made snow, enriched in 2 H (enriched 18 O water is too expensive and rare for such a field experiment), on a field site at Summit station, Greenland, a site with temperatures below 0 • C all year (at least prior to 2012).In the 4 years that followed, 2008-2011, we went back to the place four times and drilled shallow cores (from 1 to close to 4 m over the years) in which we recovered our original layer, now diffused well into the original firn surrounding it.Isotope analysis in the laboratory yielded the width of the diffused profile over the years.Together with the temperature of the layer, which was logged, and the measured density, this enabled us to compare the actually measured diffusion rate with the value from the generally used expression and parameterization by Johnsen et al. (2000).
In the following chapters, we start with the theoretical description of isotope diffusion, including our approach to numerically simulate our experiments.Next, we describe the field experiment.After that, we present the results of our measurements and discuss those extensively.We end with some conclusions and a design for a follow-up experiment.

Isotope diffusion in firn
In general, diffusion is the macroscopic description of microscopic random movements that, in combination with a gradient in the concentration of a certain constituent, cause a decrease of this gradient.The most commonly used macroscopic description originates from Fick, a 19th-century German physiologist.According to his second law, and considering only one spatial dimension, the effects of diffusion on the isotope concentration C are described as where is the diffusion coefficient, also called diffusivity, and t and z are the temporal and spatial coordinates, respectively.In our specific case, C would be the concentration of water molecules containing a 2 H isotope.In practice, however, the 2 H concentration is expressed as the deviation of this concentration from that of a reference material.This deviation is denoted by δ 2 H and is defined as where R is the abundance ratio of the rare isotope with respect to the abundant isotope: 2 H / 1 H. δ 2 H is usually expressed in per mill (‰).As the difference between concentration and ratio is very small for 2 H, to a good approximation the diffusion equation is also valid using δ 2 H. Therefore, we may change Eq. ( 1) Here m is the molar mass of water (in kg), R the gas constant (J K −1 ) and T the temperature in Kelvin.ρ f and ρ ice are the firn and ice density (kg m −3 ), respectively (ρ ice = 917 kg m −3 ).For the water vapour saturation pressure p sat we use the parameterization given by Murphy and Koop (2005): with p sat in Pa and T in Kelvin.As p sat is exponentially dependent on temperature, this parameter is the main cause for temperature dependence of the diffusion process.The other terms in Eq. ( 4), except the tortuosity τ (and m and R), are temperature-dependent as well: apart from the temperature itself, these are the ice-vapour fractionation factor α 2 and the diffusivity of deuterated water vapour in air a2 .For the most abundant water molecule 1 H 16 2 O the diffusivity in air is given in square metres per second (m 2 s −1 ) by Hall and Pruppacher (1976): where T is the temperature, T 0 is 273.15K, p is the pressure at Summit (680 hPa during summer, the time when the diffusion process is the most active) and p 0 is equal to 1013 hPa (1 atmosphere).
For water molecules containing a 2 H atom, the diffusivity is slightly lower (Merlivat, 1978) The ice-vapour fractionation factors -that is, the difference in ratio of rare and abundant isotopes in ice and vapour under equilibrium conditions -are functions of temperature and were measured by Merlivat and Nief (1967) for Deuterium: Finally, the tortuosity τ depends on the structure of the open channels in the firn.We adopt -initially -the parameterization as a function of the density of the firn that was given by Johnsen et al. (2000): This parameterization leads to increasingly high values for τ as the density ρ f approaches the density of pore close-off.
For lower densities, however, the effects due to tortuosity are assumed to be minor: according to this parameterization the value of τ varies between 1.15 and 1.25 in the density range of our experiment (300-350 kg m −3 ).Of course, this paramterization is a gross oversimplification of the real process, as it neglects the influence of varying grain sizes and shapes.Nevertheless, it seems to have served its goal reasonably well under widely varying circumstances.Diffusion decreases gradients and thus leads to an overall smoothing of the original signal.The general solution to the differential Eq. ( 3) given an initial profile δH 0 (z) is a convolution of this initial profile with a Gaussian distribution: The amount of smoothing -that is, how the values of the original profile δ 2 H 0 at positions z influence the value δ 2 H(z, t) -is determined by the width of the Gaussian curve σ 2 .The physical meaning of this width is the diffusion length, which is the average displacement of the deuterated water molecules.If the original distribution δ 2 H 0 (z, t = 0) is a Dirac distribution (infinite at z = 0, and zero everywhere else, such that its total integrated area is M), Eq. ( 10) leads to The squared value of σ is directly related to the isotopic diffusivity in firn and the elapsed time: In such an idealized case, the profile that would be recovered would show a pure Gaussian profile, and its width would be directly related to the diffusivity.The calculation of the width of such a profile would simply require the numerical integration of Eq. ( 12).For each time step, f2 needs to be calculated with the appropriate values for the variables (temperature and/or density) for that time step.
In reality, the original distribution δ 2 H 0 (z, t = 0) is of course never a Dirac function.In our experiment, however, the initial signal does resemble a Dirac function, but with a finite value for the peak value and a finite width of this peak.Furthermore, we deposit our layer of 2 H-enriched snow on a background that is not constant but that shows the natural seasonal cycle (subject to diffusion during previous years).Finally, in our experiment we sample the firn layer with a limited spatial resolution (of 1 cm).Hence we use a numerical model for the simulation of our findings, taking these complications into account (see Sect. 4.3).
Comparison of the numerically calculated σ 's with those from the field experiment enable us to test the validity of Eq. ( 4) and the terms it contains (most notably the parameterization for the tortuosity, Eq. 9).The isotope effects (Eqs.7 www.the-cryosphere.net/9/1089/2015/The Cryosphere, 9, 1089-1103, 2015 and 8) and obviously the saturation pressure of water vapour (Eq.5) are generally considered to be well known and are treated here as constants without uncertainty.Recently, Ellehoj et al. ( 2013) reinvestigated the ice-vapour fractionation factor α 2 (Eq.8) and found it to be larger, from ≈ 1 % at −15 • C to over 3 % at −40 • C.Although such changes are highly significant when studying, for example, the hydrological cycle, for our study such changes are of minor importance.

The field experiment
For the production of 2 H-enriched snow under polar field conditions, we built a snow gun installation.The snow gun itself was a small instrument designed for home garden use (CHS Snowmakers, type "Cornice").The gun produces a very fine spray of droplets which precipitate as dry, fluffy snow, provided the ambient temperature is low enough (at most −5 • C, preferably several degrees lower).We built the necessary air compressor and water pump system on a compact, gasoline-motor-driven cart.The installation was capable of producing ≈ 300 kg of snow per hour.We produced snow on an area of typically 6 m × 6 m, such that we would have ample space for drilling two-three hand cores per year for 4 years without interference.We aimed for a snow layer of 2-3 cm and allowed for loss of snow spraying outside the field, so we used about 1000 litres of water.We contained this amount of water in an inflatable children's paddling pool, which is easily transportable and also forms a good thermal isolation (we added a foam mat underneath).
The water was enriched to a level of typically δ 2 H = 1000 ‰ by adding 250 g of pure D 2 O (Sigma-Aldrich) (depending on the natural δ 2 H level of the water used).
In August 2007, we produced an area of 2 H-enriched snow in the field on pristine snow, about 2 km away from camp Summit (central Greenland, 72 • 35 N 38 • 25 W, elevation of 3216 m).The station is operated by the Americanbased CH2M HILL Polar Services (formerly Veco Polar Resources).In the summer months, there is frequent access for both people and equipment with Hercules C130 aircraft.Temperatures are always below 0 • C (at least during the years of our fieldwork).
On 8 August 2007, we produced our enriched layer in about 5 h, using ≈ 1000 L of local surface meltwater, enriched to δ 2 H = 1294 ± 3 ‰.We also dyed the water with a red food colorant, to make our produced snow layer visible.Figure 1 shows the site while producing snow.The wind speed was low that day, so most of the produced snow landed on our area (marked with poles).It was a sunny day, with temperatures reaching −5 • C, which impeded dry snow production.Therefore we produced snow at a reduced production rate.Still, the produced layer on parts of our area was ice rather than snow, especially close to the snow gun.After fin- ishing the snow production, we carefully inspected the area, such that we could try to avoid the places with ice formation when sampling in later years.
A thermistor was placed at the surface, co-located with our layer and connected to a data logger close to one of the poles.Temperatures were logged every 3 h.In this way a high-resolution continuous temperature record for our layer would be available.
Prior to our snow making, we took samples from the pristine snow layer for isotope analysis, to a depth of about 50 cm.We also performed snow density measurements, to the same depth, with 10 cm resolution.
The night after the snow was produced, the layer got covered under a few centimetres of fresh snow.
Afterwards, the depth of the snow layer was monitored by the Summit crew members every month by measuring the height of each of the five poles that marked the field; this went on until the final sampling day in 2011.At that time our snow layer was close to 3 m below the surface.
These careful snow height measurements provided us with the information needed to recover our layer in the consecutive years.In the years 2008, 2009, 2010 and 2011 (that is, 352, 643, 1102 and 1460 days after the production of the layer) we returned to Summit to drill shallow firn cores with a hand corer (Kovacs Mark II).We drilled two-three cores every year (labelled A, B, C) and made sure that we recovered the expected depth of our layer ± some 50 cm (as the depth registered at the five poles around our field scattered by 20-30 cm over the years).Figure 2 shows the depth of our layer as a function of time based on those pole height measurements, together with the points indicating the actual depth of the enriched layer (or rather the depth of the maximum δ 2 H value) as revealed by the isotope measurements later in the lab.
Still in the field, we cut the cores into 1 cm slices with a custom-built device and stored the slices in individual air- The Cryosphere, 9, 1089-1103,  tight plastic bags (Toppits Zipper).Soon after, we let the slices melt and pored them over into lockable plastic sample transport tubes (Elkay products) that had been tested for their long-term isotope integrity.
In the field, we also secured the logged temperatures of the past year, and in 2010 and 2011 we performed again 10 cm resolution density measurements, now also using our hand corer.
Back in our laboratory in Groningen, we performed δ 2 H and δ 18 O isotope ratio measurements on all samples using our routine equipment (van der Wel, 2012).The combined uncertainties were ±0.06 ‰ for δ 18 O and 0.6 to 2 ‰ for δ 2 H (depending on the level of enrichment).In all of the total of 10 cores drilled over the years, we found our enriched layer, close to the depth expected based on the pole height measurements (Fig. 2).

Density and temperature
For the simulation of the diffusion of our enriched layer, reliable values for both the temperature and the density are the most important input values.Figure 3 shows the density measurements that we performed over the years, all measured close to the area of the enriched layer, grouped in a single plot.The depths have been shifted (using the information shown in Fig. 2) such that our enriched layer is at depth zero.The data show that our enriched layer, deposited at the end of summer, is on top of a layer with lower density than the preceding and following winters.This summer-winter effect is beautifully demonstrated by Albert and Shultz (2002) from Summit in 2000, and our data are in agreement with their findings (shown in their Fig. 2).Based on their and our data we use an initial density of 300 kg m −3 for our diffusion calculations, with a gradual increase of 10 (kg m −3 ) yr −1 .
The temperature registration of the thermocouple, at the same depth as the layer, is shown in Fig. 4. Unfortunately, in spite of our efforts, two larger parts of the total temperature profile were lost.Figure 4 shows the interpolations that we made.We estimate the extra uncertainty in the diffusion calculations due to this omission to be minor.Fortunately, the first full year of data has been recovered.This is the part when the layer is still so shallow that the diurnal temperature cycle (which we capture by our 3-hourly temperature sampling) is still noticeable (see insert in the figure for the first month).As the diffusion rate is exponentially dependent on temperature, capturing this first period in detail is crucial for the results of the numerical simulations.

Field diffusion profiles
For each sampling year, two-three records for both δ 2 H and δ 18 O were measured (labelled A, B, C). δ 2 H contains the crucial diffusion information: the broadened (and weakened) profile around our original layer of enriched δ 2 H.The quality of our collected profiles was variable.Some of the profiles showed one or two samples (= cm) that had δ 2 H values almost as high as the original enriched water, whereas all other samples were close to the natural values.We attribute this to ice formation during the snow production, reducing diffusion rates dramatically.Fortunately, for every year we also had profiles without such irregularities, which showed a clear, Gaussian-shaped profile above background.Table 1.The results for the diffusion length σ and for the net peak height for all profiles.The increase of σ as a function of time is clearly visible.The uncertainties in the second column are those from the fitting procedure.The final combined uncertainties in the results are presented in the third column.Except for 2008, all measured values per year agree within their uncertainties.The net peak height is also dependent on the initial thickness of the enriched layer.Therefore the found peak height is variable within and between years.Figure 5 shows two of the δ 2 H profiles, 2008B, and 2011A, respectively.The effect of diffusion is directly visible, both in the width of the peak, and in its height.For the quantitative determination of both, however, we need to correct for the natural δ 2 H seasonal cycle that interferes with the diffused pattern of our original enriched layer.We used the δ 18 O profile to reconstruct the natural δ 2 H seasonal cycle.δ 2 H and δ 18 O in precipitation show both a very similar seasonal cycle, with the amplitude of the δ 2 H cycle being around 8 times as large as that of δ 18 O.Contrary to that of δ 2 H, the δ 18 O seasonal cycle is not influenced by our layer: the water used was in fact recent snow at Summit, with δ 18 O ≈ −30 ‰ very close to the value of the top layer of our field.In Fig. 5, δ 18 O is shown as well, with scale ratio 1 : 8 with respect to the δ 2 H scale.For the reconstruction of the natural δ 2 H signal the δ 2 H − δ 18 O ratio for all our 10 profiles was fitted individually, by using the flanks of the profiles.Subsequently, we corrected our measurements for this reconstruction of the natural δ 2 H signal, thereby obtaining the net diffusion profile.
For each of the 4 years, we had two profiles available (and even 3 for 2008 and 2010).Only one of these 10 profiles (2008A) was not useful: the deposited layer consisted only of ice and diffusion had hardly taken place.Figure 6 shows all other net δ 2 H profiles, together with the Gaussian fits, after subtracting the background signal.The width of the fit, which is the diffusion length σ (see Eqs. 10 and 11), is also indicated.
Not all profiles are of equal quality: half of them showed the presence of ice inside our deposited layer, visible through one isolated high δ 2 H value in the profile (on most occasions the ice layers had already been noticed in the field); as those points are not representative for the diffusion process, we discarded them.This happened for profiles 2008B, 2009B, 2010B, 2010C and 2011A.Furthermore, some points had to be discarded that resulted from contamination of samples with snow/firn from other depths that happened during the coring process (such contamination was also visible in the δ 18 O signal).Discarded points are shown in the plots in brown.Table 1 shows the results for the diffusion length σ and the net peak height for all profiles.The increase of σ as a function of time is clearly visible.Contrary to σ , the net peak height is not only dependent on the diffusion time, but also on the initial thickness of the enriched layer.Therefore the found peak height is expected to be variable within and between years.
The uncertainties given in the second column in Table 1 are those from the fitting procedure.While they give a good  indication for the fit quality, the final combined uncertainty in the results is, of course, considerably higher.The main experimental uncertainty lies in the representation of the "z axis", the depth.We estimate this error to be ±3 % of the value, leading to an error in σ of about 0.10.The icy character of our deposited layer in some profiles form another principal source of error: although on both sides of such an ice layer the firn diffusion process takes place, and we can thus use those profiles for a σ measurement, the width of the fitted Gaussian curve will still be influenced by the presence of the icy character of the original layer itself.Therefore, we have increased the uncertainties for such profiles to ±0.25 cm.The www.the-cryosphere.net/9/1089/2015/The Cryosphere, 9, 1089-1103, 2015 uncertainty caused by the δ 18 O-based background correction is negligible.The final attributed errors are given in the "uncertainty" column.Except for 2008, all measured values per year agree within these uncertainties.As each year had at least one core with, and at least one core without, ice in our deposited layer, the fact that their diffusion lengths agree with each other shows that these ice layers did not influence the diffusion length significantly in this experiment.

Comparison with the numerical simulation
The simplest way of simulating our experiment is to numerically integrate Eq. ( 12) using the known temperature and density as a function of time (Eq.4).However, the real experimental situation is more complicated.To accurately simulate the experimental situation, we first calculated the diffused δ 2 H pattern as a function of time from the original δ 2 H 0 (z, t = 0) pattern around our enriched layer with an added "pulse" of enriched δ 2 H.We know the value of this enriched δ 2 H (1294 ‰), but the thickness of the layer is unknown and variable.Therefore we calculated the profiles for three initial layer thicknesses: 6, 18 and 30 mm.As the next step, we corrected the diffused pattern for the slight compaction that occurred (inversely proportional to the small increase in density), and we sampled the diffused patterns with the spatial resolution of the experiment (1 cm).Then, we simulated the correction for the natural δ 2 H seasonal cycle using the also diffused and sampled δ 18 O profile, in the same way They are compared to the numerical simulations with initial layer thicknesses of 6, 18 and 30 mm, for which the diffusion length is fitted to the experimental points in Fig. 7.All experimental points are in the range spanned by the numerical calculations.The values for the ice layers that we removed from our profiles are also indicated (half-filled squares): all but one lie outside the possible range for diffused firn profiles, identifying them once more as ice layers.
as we did with the experimental profiles.Finally, we fitted the net δ 2 H profile with a Gaussian function.
Figure 7 shows the results for σ achieved this way, as well as the σ from the direct integration of Eq. ( 12).The differences between the numerical calculations at variable initial thickness are quite small, especially for the values of σ for later years, indicating that the effects of sampling and the background correction are minimal.Figure 7 also contains the experimental diffusion lengths and thus embodies the main results of this work.Clearly, there is a systematic mismatch between the experimental results and the numerical simulations, increasing with time.To fit the data, the simulated curves need to be ≈ 25 % lower in the first year, up to ≈ 40 % in the final year.This implies a lowering of f2 up to a factor of 2.5 (as σ is proportional to the square root of f2 ).This lowered fit curve is also shown in Fig. 7. (The black dotted line in between will be described in the discussion section.)All in all, Fig. 7 suggests that either there is an experimental flaw or else one or more parameters included in Eq. ( 4), or Eq. ( 4) as a whole, are not adequate.Below we discuss various possibilities influencing the average value of σ and its dependence on time.
The full numerical procedure also results in a peak height, which at any point of time is proportional to the width of the initial enriched δ 2 H pulse. Indicated in Fig. 8 (red dots) are the actual peak height fits of the profiles (see Table 1).They are compared to the numerical simulations for layers with initial thicknesses 6, 18 and 30 mm, for which the diffusion length is fitted to the experimental points in Fig. 7 (the red dotted line).All experimental points are in the range spanned by the numerical calculations.The values for the ice layers that we removed from our profiles are also indicated: The Cryosphere, 9, 1089-1103, 2015 www.the-cryosphere.net/9/1089/2015/ the position of four out of five of them in this plot clearly corroborates them as ice layers, as an initial layer thickness of 30 mm is about the maximum realistic value for our snow layer.

Discussion
The considerable discrepancy between our experimental results and the numerical calculations based on Johnsen et al. (2000) came as a surprise.The theoretical description of the firn diffusion process by Johnsen et al. (2000), a further development of work by Johnsen (1977) and Whillans and Grootes (1985), has been used for the description (and back correction) of diffusion in many ice core projects.
The large majority of the papers describes, or "backcorrects", the influence of firn diffusion as it is recorded in the ice below pore close-off.There, the ice carries the result of firn diffusion integrated over the entire firn phase.Examples of such work, restricted to Greenland, are Vinther et al. (2006Vinther et al. ( , 2010)), Masson-Delmotte et al. (2005), Andersen et al. (2006b), Jouzel et al. (1997), White et al. (1997) and Simonsen et al. (2011).The last publication concentrates on the so-called differential diffusion, the difference in diffusion between δ 18 O and δ 2 H, which is only dependent on the temperature of the firn while diffusion takes place.This idea was in fact the main subject of Johnsen et al. (2000).Below we will dicuss three possible causes for the discrepancy.They are (1) our experimental conditions, especially the formed ice in the deposited layer; (2) a considerable influence of tortuosity already in the uppermost metres of the firn, contrary to the assumptions in Johnsen et al. (2000); and (3) invalidity of the assumption that no gradient in isotopic composition builds up within the firn grains.

Experimental conditions
The occurrence of an ice layer can practically block the diffusion process.We do, however, firmly believe that our results have not been seriously influenced by ice formation.Each year contained both a profile with and one without the indication of ice formation inside our deposited layer.Yet, for each of the 4 years that we sampled, the results for the diffusion length agree very well.As the occurrence of an ice layer almost stops the diffusion process (see e.g.van der Wel et al., 2011b), one would expect large scattering within and between years if ice formation inside our deposited layer indeed played an important role.The fact that they do not can be explained by the fact that the water vapour from such an ice layer will immediately encounter natural snow layers in which the diffusion process occurs naturally.Only the ice layer itself will continue to contain a high level of enrichment and in the end produces a δ 2 H value that needs to be excluded from the fit to the data, as we did.
Even in the absence of ice, the density of our artificial snow layer is probably higher than that of fresh Summit snow.When trying to fit the results of Fig. 7 using higher densities we find that we would need a more or less plausible density of around 380 kg m −3 for the first year, but increasingly higher values for the subsequent years, up to 520 kg m −3 for the 2011 results.Such compaction of a layer, initially already denser than its surroundings, in just 4 years is unrealistic.Furthermore, again the diffusion process will immediately encounter natural snow as soon as the process starts.So, in the course of the years, with diffusion lengths getting larger and the signal profile getting dominated by the region outside the original layer, one can expect any initial effect of higher density to weaken.However, we observe the opposite: the deviation between our experimental results and the simulations based on Johnsen et al. (2000) increases in the course of the years.

Tortuosity
The diffusion length σ is inversely proportional to the square root of the tortuosity τ .If the discrepancy between our results and the numerical simulation were entirely due to higher tortuosity in our experiment than the range of 1.15-1.3given by Eq. ( 9), we would need the tortuosity to be between 2.5 and 3. To see to what extent this would be plausible, we have gathered relevant information from the literature describing both firn isotope diffusion and gas diffusion.To facilitate a proper comparison, we first have to define the tortuosity in an unambiguous manner: Here a is the diffusivity of the compound (water vapour in our case or, more precisely, deuterated water vapour) through a certain area of free air, and f the effective diffusivity through that same area, but now filled with firn.The porosity ϕ accounts for the effective open area available for the diffusion process, whereas the tortuosity accounts for the shape of the air channels.In the case of perfectly straight air channels, τ would be 1.
The Whillans and Grootes model (Whillans and Grootes, 1985) was the first to describe firn diffusion in a detailed manner, but they did not include the influence of tortuosity.Instead of the porosity, they included the density at pore close-off ρ c : Using this ϕ * is equivalent to using Eq. ( 13), with As the difference between ρ ice and ρ c is rather small (Whillans and Grootes (1985)  fective tortuosity values remain close to 1 except for densities approaching ρ c .Cuffey and Steig (1998) performed a detailed -that is, time-resolved -study of firn diffusion showing the dampening of the seasonal cycle in δ 18 O in the shallow GISP-B core (based on the data published by Stuiver and Grootes (2000) and Stuiver et al. (1995)) and compared that to the Whillans and Grootes (1985) model (taking the low atmospheric pressure at Summit into account).They concluded that the model agreed well for the upper metres (starting, however, from a depth of 1.5 m) but that the diffusion effects produced by the model were too large for larger depths.The diffusivity of the model needed to be decreased by a factor of about 1.7 to match the data.They built this into the model by decreasing the maximum density at which firn diffusion still occurs, from 830 down to (a fitted best value of) 730 kg m −3 .Using Eq. ( 15) their results can again be expressed as using a tortuosity factor that is now considerably higher than in the original Whillans and Grootes (1985) model.However, even in this study the depth resolution is limited to the length of one seasonal cycle (typically 50 cm); they ignored the top 1.5 m, and their numerical procedure concentrated on the amplitude of the seasonal cycle only, not taking into account for example the diffusion differences between original summer and winter snow.The temperature driving the diffusion process is parameterized.Johnsen et al. (2000) modified the Whillans and Grootes (1985) model, the main difference being the explicit introduction of the tortuosity factor.Other differences, with relatively small (< 5 %) influence, are a different parameterization of the water vapour diffusivity through free air, and the more complete treatment of the isotope effects.The tortuosity that Johnsen et al. (2000) used (Eq.9) is based on a fit to gas diffusion measurements by Schwander et al. (1988), performed on firn samples from Siple station, Antarctica.The density range of these measurements was between 500 and 750 kg m 3 .Schwander et al. (1988) present the tortuosity as it is defined in Eq. ( 13) in their Fig. 5. Tortuosity values they found increased with the density from 2 to 7. They also provided a table (their Table 1) that may give rise to some confusion.The effective diffusivity that is given there is in fact f /φ ("the flux per unit cross section in the open pores") (J.Schwander, personal communication, 2014).So, dividing the diffusivity given in that table by the open air diffusivity a directly yields 1/τ .Their results for tortuosity were generally, although coarsely, confirmed by Jean-Baptiste et al. (1998), who did the first diffusion measurements on deuterium isotopes in firn: they measured the isotope diffusion around the connection of firns with distinctly different isotopic composition.The firn was, in fact, crushed ice.They used densities between 580 and 760 kg m −3 .Both the Schwander et al. (1988) results and those by Jean-Baptiste et al. ( 1998) have thus been performed for higher densities only.Using the Schwander et al. (1988) parameterization by Johnsen et al. (2000) for our experiment, with densities varying from 300 to 350 kg m 3 , means a substantial extrapolation.
More recent measurements of firn diffusivity are presented in studies by Pohjola et al. (2007) and by van der Wel et al. (2011a).Both studies have improved on the work by Jean-Baptiste et al. (1998), since they have identified the interface between the two stacks with different isotopic composition as the weak spot of such experiments.By connecting several layers of different thickness, they could identify possible problems with these interfaces, for example when such an interface was much more porous than the bulk material.Due to these effects, the results by Pohjola et al. (2007), at a density range of 480-500 kg m −3 , did not produce consistent results for the tortuosity.Van der Wel et al. (2011a), however, managed to make an ideal multi-layer snow stack (produced with the same snow gun that we used in the present study).At a density of 415 kg m −3 , they were able to fit their diffused isotope profiles using the expression by Johnsen et al. (2000), thereby finding a tortuosity of 1.18 ± 0.08 (compared to the value 1.36 that follows from Eq. 9).
Whereas this latter piece of information indicates tortuosity values hardly above 1 for low densities, various gas diffusion measurements show considerably higher values.Fabre et al. (2000) performed gas diffusivity measurements on site in Vostok, Antarctica, and at an alpine site (Col du Dome).They express their results in the same fashion as Schwander et al. (1988) (they also show their results for comparison) and also give various model results for the tortuosity.Albert and Shultz (2002) performed detailed gas diffusivity and permeability measurements on the top 10 m of snow and firn at Summit station.They mention the tortuosity, defined in their paper as the reciprocal from our Eq.( 13) as a side result, and quote the value of ∼ 0.5 for the top layer of the snow.When looking at their data, however, it seems that they used a different definition for the tortuosity and actually determined the ratio f / a to be 0.5.In more recent work by the same group, they avoid the term tortuosity altogether and instead report on the f / a ratio.Adolph and Albert (2013) describe an improved way to measure gas diffusivity through firn, and they report a series of diffusion measurements performed on firn from Summit.From the results in their Table 1 (in which they also quote their previous result from 2002) we can deduce the tortuosity as defined in Eq. ( 13).A subsequent paper by the same authors (Adolph and Albert, 2014) reports on even more diffusivity measurements (and includes the ones reported in 2013).
Recently several large firn gas transport studies were published, one of them (Buizert et al., 2012) concentrating on the NEEM site in northern Greenland.These authors tuned six firn air transport models to firn concentration measurements of a set of 10 reference trace gases.Whereas the fit quality of the tracer concentrations is high, for the upper 4.5 m the fit is underdetermined, and the spread of the molecular diffusivity profiles for CO 2 is large.Furthermore, convective mixing The Cryosphere, 9, 1089-1103, 2015 www.the-cryosphere.net/9/1089/2015/plays a role in modelling these upper 4.5 m, a process that influences gas transport far more than the firn itself (see below).
In direct firn gas diffusion experiments convective mixing is avoided, making such results more useful for describing firn diffusivity.
We give an overview of the tortuosity from all these results, according to our Eq.( 13), in Fig. 9.This figure includes the tortuosity by Cuffey and Steig (1998), which follows from Eq. ( 15), and the parameterization for the Schwander et al. (1988) data used by Johnsen et al. (2000) (and thus also by us in the previous chapters).Furthermore, as a lower limit, we have included the theoretical result by Weissberg (1963), which he derived for spheres that can partly overlap or even fuse: Using Eq. ( 13) this is equivalent to Looking at Fig. 9, we observe large scatter indeed.In the higher-density region the tendency towards higher values for τ is clear, but there is a considerable discrepancy between the Schwander et al.However, the Adolph and Albert ( 2014) results also suggest higher values for τ at lower densities.This is in contrast with the only firn vapour measurement at low densities by van der Wel et al. (2011a) that gives a value that is lower than the Weissberg (1963) theoretical model.In contrast to that, the Adolph and Albert (2014) results suggest that a higher tortuosity in the low-density region of ≈ 1.5 is probably a better choice than the low values given by the extrapolation of the Johnsen et al. (2000) parametrization.The 1 Figure 6 in Adolph and Albert (2013) suggests the opposite.This is, however, because the authors interpreted the results in Table 1 of Schwander et al. (1988) as f , whereas they are in fact f /ϕ (A. C. Adolph, and J. Schwander, personal communication, 2014).13).The symbols are the results of the various measurements: the red squares are laboratory firn diffusion measurements (by Jean-Baptiste et al., 1998, andvan der Wel et al., 2011a); the blue circles CO 2 and O 2 gas diffusion measurements through firn from Siple dome, Antarctica (Schwander et al., 1988); and the green triangles SF 6 gas diffusion measurements through firn from Summit (Adolph andAlbert, 2013, 2014).Finally pale brown triangles are the results by Fabre et al. (2000) for SF 6 gas diffusion in Vostok and Col du Dome.The black dashdotted line is the parameterization for the Schwander et al. (1988) data in the Johnsen et al. (2000) model.The grey dashed curve is the tortuosity that is equivalent to the fit that Cuffey and Steig (1998) made to isotope seasonal cycles in the Summit firn layer.Furthermore, as a lower limit, we have included the theoretical result, interpreted as τ , by Weissberg (1963), which he derived for spheres that can partly overlap or even fuse.For discussion of the results, see text.
value range of 2.5-3, however, that we need to fit our data, is not supported by any data in the low-density range.In Fig. 7, we show the numerical results for σ using a fixed tortuosity of 1.6.Whereas for the first year agreement is reasonable, in the following years the deviation increases: in the experiments, diffusion is slowing down, and this is unlikely to be caused by tortuosity (or density) increases.
The clearest conclusion of all, however, is that the parameterization of τ as a function of density (or porosity) is an oversimplification.Albert and Shultz (2002) show the structural changes of the fresh snow in its first years: whereas the density hardly changes, grain size rapidly grows, and the permeability (and likely also the diffusivity) increases.In that paper, they report a single diffusivity result (included in Fig. 9) on the top 20 cm of the firn, yielding τ = 1.36 for ρ = 326 kg m −3 .This suggests that, for the youngest firn, time since deposit is a more important parameter than density.Furthermore, τ will have a different course in time for winter than for summer snow.All in all, we conclude that for the simulation of our experiment, choosing a τ of ≈ 1.6 would be a fair choice given the data available, but a value range of τ = 2.5-3, needed to fit our data using the Johnsen et al. (2000) simulation, is not plausible.

Isotope homogeneity within the firn grains
In the model by Johnsen et al. (2000), several assumptions have been made: -the effects of firn ventilation are negligible, -there is continuous isotopic equilibrium between the ice grain surfaces and the vapour, -the ice grains themselves remain isotopically homogeneous.
One or more of these issues have been addressed by several authors, among whom Whillans and Grootes (1985), Jean-Baptiste et al. ( 1998), Johnsen et al. (2000) themselves and Neumann and Waddington (2004).The latter paper describes a very detailed numerical model in which in the first place the influence of firn ventilation is quantified.They conclude that isotope exchange in the upper few metres is more rapid than follows from models such as that of Whillans and Grootes (1985) and Johnsen et al. (2000).However, they also state that the ventilation process is especially important in low-accumulation zones, such as the Antarctic For the Summit site, its influence will probably only be marginal.Moreover, whereas firn ventilation might lead to changes in the overall isotopic composition, its character will not be the same as diffusion; especially it will not influence the diffusion pattern, and thus the diffusion length fit, of our enriched layer.
The model by Neumann and Waddington (2004) also allows for a disequilibrium between the ice grain surface and the vapour.Although the isotope exchange rate is not well known, they conclude that the ice phase is not in isotopic equilibrium with the vapour at any depth in the firn.This effect would slow down the influence of diffusion.Somewhat surprisingly, however, these authors assume the grains themselves to be and remain isotopically homogenous.This is in fact the point that the other authors touch upon.Isotope diffusion in ice, and thus inside the grains, is 10 orders of magnitude smaller than vapour diffusion through air.At −20 • C, for example, the solid-ice diffusivity is about 1 × 10 −15 m 2 s −1 (Whillans and Grootes, 1985), whereas vapour diffusion through air, according to Eq. ( 6) (Hall and Pruppacher, 1976), yields 2.7 × 10 −5 m 2 s −1 at Summit.In contrast, the water molecules spend 5 to 6 orders of magnitude more time in the solid than in the vapour phase (depending both on temperature and density of the firn).Both Whillans and Grootes (1985) and Johnsen et al. (2000) have incorporated this into their model by dividing the free vapour diffusion rate by this residence time ratio.Nevertheless, the solid-ice diffusivity remains some 5 orders of magnitude slower than the effective vapour diffusion.Whillans and Grootes (1985) have investigated this and concluded that, given the average size of the firn grains, the isotopic homogeneity assumption is valid.Jean-Baptiste et al. (1998), however, show in a numerical model that their model description of their firn diffusion experiment would indeed be influenced for grain sizes of 1mm and larger.Johnsen et al. (2000) conclude that grain homogeneity will not occur based on ice diffusion alone "for the coarse grained (2 mm) summer layers".As they on the other hand conclude from observations that the isotopic seasonal cycle can disappear completely due to diffusion, they propose grain boundary migration as a different mechanism for more rapid grain isotope homogenization.
The ratio between the effective firn vapour diffusion and solid-ice diffusion is the largest for low densities, as then the solid-to-vapour residence time ratio is the lowest.These are the circumstances of our experiment.Furthermore, the detailed microstructure experiments by Albert and Shultz (2002) show substantial growth of grain size in the first years after deposition.Together, these circumstances will probably cause a substantial inhomogeneity in the ice grains, thereby slowing down diffusion to below the rate described by Eq. ( 4).

Conclusions and outlook
With this experiment, we have observed isotope diffusion in the natural setting of Summit in the first 4 years after deposit down to about 3 m depth.The idea of determining the diffusion length as the Gaussian width of an initial thin layer substantially enriched in deuterated water worked very well indeed.The results, however, indicate a substantially lower diffusivity than expected based on the well-established model by Johnsen et al. (2000).Our attempt to explain this difference brought us to the following conclusions: -Although we can not be fully sure whether the characteristics of the enriched layer itself, with some local ice formation, has slowed down the diffusivity, it is very likely that the diffusion lengths we have obtained resemble the true diffusivity of deuterated water in the upper layers of Summit firn, because the diffusion takes place in the original snow layers after some time.
-Tortuosity is in general poorly characterized.Several of the firn and gas diffusion experiments over the years lead to a very scattered total picture.The parametrization used by Johnsen et al. (2000) is probably not correct for Summit; based on recent gas diffusion experiments by Adolph and Albert (2014), tortuosity is probably considerably larger in the uppermost layers, but in contrast not as large in the deeper firn.Density is a poor measure for tortuosity, certainly in the upper metres of The Cryosphere, 9, 1089-1103, 2015 www.the-cryosphere.net/9/1089/2015/firn, and considerable differences between summer and winter layers are likely to exist.
-Nevertheless, the discrepancy between our results and the Johnsen et al. (2000) model cannot be explained by higher tortuosity alone, as the value that would be needed to fit our data is definitely outside the plausible range.
-It is likely that isotopic inhomogeneity exists within the ice grains in the firn, as the vapour diffusion process is orders of magnitude faster than solid-ice diffusion.This effect is only partly compensated by the much longer residence time of the molecules in the solid phase.This inhomogeneity, and its slowing effect on diffusion, depends critically on grain size.In the first years after snow deposition, grains tend to grow (Albert and Shultz, 2002), thereby effectively slowing down diffusion.Of the three possible causes for the discrepancy between our data and the simulations, isotopic inhomogeneity is thus the most plausible: it would explain why gas diffusion measurements (and thus also the parameterization used by Johnsen et al., 2000) are not entirely valid for firn vapour diffusion.In the first stage of the diffusion process, the part that our experiment monitors, ice grain inhomogeneity would slow down the diffusion process, but in a later stage the inhomogeneity would diminish or even disappear again.Combined with lower values for the tortuosity at greater density (such as for instance those by Adolph and Albert, 2014)  -Although the model by Johnsen et al. (2000) has been used extensively and successfully, this has to our knowledge never been done in the amount of detail we present.Rather, it has been used to describe the integrated diffusion over the entire firn phase, expressed as the total diffusion length caused by it.This integrated diffusion length is a rather forgiving parameter: it does not matter whether the diffusion length grows rapidly initially and then slows down, or if it increases more steadily over the years.Moreover, due to compaction the diffusion length even decreases again, although the diffusion effects (such as the decrease of the amplitude of the seasonal cycle) of course are not lessened.According to the model, the major part of the diffusion length is built up in just the upper 10 m (see also Simonsen et al., 2011), and the maximum values are achieved around 30 m depth; from then onwards, the compaction leads to a gradual decrease in diffusion length.These features of the model have, to the best of our knowledge, never been checked experimentally.
An experiment to perform such a decisive test would consist of two parts: first is the high-resolution (typically 2 cm) isotope measurement of the upper ≈ 30 m of a firn core.With a modern, laser-based isotope measurement system on-site, this would probably be feasible in just one summer season.If designed carefully, this set-up delivers the density of the firn as well.Second is reconstruction of the input function: the temperature, the precipitation events and -ideally -their isotopic composition.In this way, the "virtual ice core" approach (van der Wel et al., 2011b) can be followed to reconstruct the non-diffused firn profile, and, by comparing that profile to the data, diffusivity can be calculated with high temporal/spatial resolution.The results can then be compared to the output of the Johnsen et al. (2000) model, and to a model containing both the firn diffusion and the diffusion inside the grains.A valuable additional measurement would be isotope measurements on vapour in firn air at various depths.In this way the isotopic (dis)equilibrium could directly be assessed.
Summit would be the ideal spot for such an experiment: it has been in operation since 1988, and surface temperature and precipitation amounts have been logged since then.Furthermore, many scientific experiments run at Summit each year, many of them including stable isotope measurements.
Such a detailed study would finally enable us to describe isotope diffusion in firn in a reliable, quantitative way.The results of such a study would especially be crucial for the use of differential diffusion to reconstruct palaeotemperatures: the way the diffusion process behaves through the firn layer determines the weighted average of the temperature that is conserved in the differential diffusion signal.If the upperfew-metre diffusion is less prominent than the Johnsen et al. (2000) model suggests, this weighted average would be far less sensitive to the high summer day temperatures than is assumed at present.

Figure 1 .
Figure 1.The production of the 2 H-labelled layer of snow at Summit, 8 August 2007.From left to right the air compressor/water pump rack, the inflatable water container and the snow gun.

Figure 2 .
Figure 2. The depth of the enriched snow layer as a function of time, based on measuring the height of each of the five poles that mark the field (lines), together with the points indicating the actual depth of the enriched layer as revealed by the isotope measurements later in the lab (two to three points per year).

Figure 3 .
Figure 3. Density measurements, with depth resolution of 10 cm, performed at our site in the years 2007, 2010 and 2011.The depths have been shifted such that our enriched layer (with an estimated thickness of typically 2 cm) is at relative depth zero.Our enriched layer, deposited at the end of the summer of 2007, is on top of a summer layer with lower density than the preceding and following winter layers.

Figure 4 .
Figure 4.The temperature registration of the thermocouple, at the same depth as the layer.The dashed lines are interpolations for the times we did not have temperature measurements available.The insert shows the first month of data in detail.

Figure 5 .
Figure 5. Two examples of δ 2 H profiles obtained in this work, from 2008 (core B) and from 2011 (core A).The effect of diffusion is directly visible, both in the width of the peak and in its height.The δ 18 O signal is shown in blue, with scale ratio 1 : 8 with respect to the δ 2 H scale.

Figure 6 .
Figure 6.The nine net δ 2 H profiles for the four consecutive years, together with the Gaussian fits.The width of the fit, which is the diffusion length σ , is listed in the plots.The brown circles are measurements discarded from the fit for various reasons (see text).The increase of σ over the years is clearly visible.

Figure 7 .
Figure 7.The results for the diffusion length σ .The points are the experimental results.The set of four higher-lying curves (one black, three in shades of blue) is the result of the numerical simulation using the Johnsen et al. (2000) model.The black one is the direct calculation of σ ; the three blue ones result from the full numerical procedure starting with an enriched layer with 6, 18 and 30 mm initial thickness.Clearly, there is a systematic mismatch between the experimental results and the numerical simulations, increasing with time.The fit curve through the measurements is constructed by lowering the diffusivity by 25 % in the first year, up to 40 % in the last.Finally, the black dotted line shows the numerical result using a fixed, higher value for τ (1.6).

Figure 8 .
Figure 8.The actual peak height fits of the profiles (red circles).They are compared to the numerical simulations with initial layer thicknesses of 6, 18 and 30 mm, for which the diffusion length is fitted to the experimental points in Fig.7.All experimental points are in the range spanned by the numerical calculations.The values for the ice layers that we removed from our profiles are also indicated (half-filled squares): all but one lie outside the possible range for diffused firn profiles, identifying them once more as ice layers.
used 830 kg m −3 for ρ c ), ef- (1988) values (based on CO 2 and O 2 diffusion through firn from Siple Dome, Antarctica) and the Fabre et al. (2000) ones (Vostok, Antarctica, and the Col du Dome alpine site based on SF 6 and CO 2 diffusion) on the one hand and the results by Adolph and Albert (2014) (using SF 6 diffusion through firn from Summit station) on the other.Although the Adolph and Albert (2014) results show some higher values, in general their results for τ are much lower 1 , and thus the diffusion process would occur more rapidly.The real firn vapour diffusivity experiments by Jean-Baptiste et al. (1998), especially the highest density measurement, seem to corroborate the Schwander et al. (1988) and Fabreet al. (2000)  results.Furthermore, as is clear from the work ofCuffey and Steig (1998), the total firn diffusion process on Summit can be described using their parameterization, which is not very different from theSchwander et al. (1988) measurements, and theJohnsen et al. (2000) parameterization.

Figure 9 .
Figure9.Results of various tortuosity measurements from the literature, with τ defined according to our Eq.(13).The symbols are the results of the various measurements: the red squares are laboratory firn diffusion measurements (byJean-Baptiste et al., 1998, andvan der Wel et al., 2011a); the blue circles CO 2 and O 2 gas diffusion measurements through firn from Siple dome, Antarctica(Schwander et al., 1988); and the green triangles SF 6 gas diffusion measurements through firn from Summit(Adolph and Albert, 2013, 2014).Finally pale brown triangles are the results byFabre et al. (2000) for SF 6 gas diffusion in Vostok and Col du Dome.The black dashdotted line is the parameterization for theSchwander et al. (1988) data in theJohnsen et al. (2000) model.The grey dashed curve is the tortuosity that is equivalent to the fit thatCuffey and Steig (1998) made to isotope seasonal cycles in the Summit firn layer.Furthermore, as a lower limit, we have included the theoretical result, interpreted as τ , byWeissberg (1963), which he derived for spheres that can partly overlap or even fuse.For discussion of the results, see text.
the total integrated diffusivity would still fit to the well-known isotope signals in the ice.To prove (or disprove) this hypothesis, a new model framework needs to be developed to incorporate grain inhomogeneity into the model.Jean-Baptiste et al. (1998) and especially Neumann and Waddington (2004) have shown pathways towards such models.