Dispersion in deep polar firn driven by synoptic-scale surface pressure variability

. Commonly, three mechanisms of ﬁrn air transport are distinguished: molecular diffusion, advection, and near-surface convective mixing. Here we identify and describe a fourth mechanism, namely dispersion driven by synoptic-scale surface pressure variability (or barometric pumping). We use published gas chromatography experiments on ﬁrn samples to derive the along-ﬂow dispersivity of ﬁrn, and combine this dispersivity with a dynamical air pressure propagation model forced by surface air pressure time series to estimate the magnitude of dispersive mixing in the ﬁrn. We show that dispersion dominates mixing within the ﬁrn lock-in zone. Trace gas concentrations measured in ﬁrn air samples from various polar sites conﬁrm that dispersive mixing oc-curs. Including dispersive mixing in a ﬁrn air transport model suggests that our theoretical estimates have the correct order of magnitude, yet may overestimate the true dispersion. We further show that strong barometric pumping, such as at the Law Dome site, may reduce the gravitational enrichment of δ 15 N–N 2 and other tracers below gravitational equilibrium, questioning the traditional deﬁnition of the lock-in depth as the depth where δ 15 N enrichment ceases. Last, we propose that 86 Kr excess may


Introduction
The firn layer is the upper 50-120 m of consolidating snow found in the accumulation zones of ice sheets and glaciers.Within this perennial snowpack a network of connected pores exists that facilitates the movement of air.The firn layer is a mixed blessing.On the one hand it complicates the in-terpretation of ice core records via the gas age-ice age difference (or age; Schwander and Stauffer, 1984), waterisotope diffusion (Johnsen, 1977), the broadening of the gas age distribution (Schwander et al., 1993), diffusive isotopic fractionation of trace gases (Trudinger et al., 1997;Buizert et al., 2013), and non-atmospheric gas variations originating from layered bubble trapping (Etheridge et al., 1992;Rhodes et al., 2016).On the other hand, the firn provides a valuable archive of old air (Etheridge et al., 1996;Battle et al., 1996), and the characteristics of firn air movement give rise to additional signals that can be used as proxies for local temperature change (Severinghaus et al., 1998;Leuenberger et al., 1999) and surface elevation (Raynaud and Lorius, 1973), or as tools for ice core dating (Bender, 2002;Kawamura et al., 2007).It is therefore clear that a complete scientific understanding of firn air transport is critical to both the correct interpretation of the ice core record, and for utilizing the unique possibilities offered by the firn-derived proxies.
Commonly, firn air transport models include three mechanisms of air movement (Trudinger et al., 1997;Rommelaere et al., 1997;Buizert et al., 2012).The first is downward advection with the ice matrix, which operates along the full length of the firn column.In a Eulerian frame of reference, the downward velocity of the air (w air ) is smaller than that of the ice itself (w ice ) due to a back flux of air from pore compression (Buizert, 2013).The second mechanism is molecular diffusion in the open pores.As the tortuosity of the pore space increases with depth, the effective diffusivity decreases; at the so-called lock-in depth, molecular diffusion effectively ceases.Molecular diffusion in the vertical direction leads to gravitational enrichment in δ 15 N-N 2 and other tracers.The third mechanism is convective mixing, an Published by Copernicus Publications on behalf of the European Geosciences Union.
C. Buizert and J. P. Severinghaus: Dispersion in firn by pressure variability umbrella term for several phenomena that vigorously ventilate the upper few meters of the firn column (Powers et al., 1985;Colbeck, 1989;Kawamura et al., 2006;Severinghaus et al., 2010).Convective mixing is commonly implemented in models as an eddy-type diffusivity that is equal in magnitude for all gases, and does not enrich gravitationally (Kawamura et al., 2006).
Diffusion effectively ceases at the lock-in depth, and consequently there is no gravitational enrichment within the lock-in zone (Sowers et al., 1992).However, detailed studies at the NEEM site in northern Greenland and the Megadunes site in Antarctica suggest continued vertical mixing within the lock-in zone (Buizert et al., 2012;Severinghaus et al., 2010).Figure 1 shows the 14 CO 2 concentrations in NEEM firn air samples (black dots with error bar), which reflect the atmospheric "bomb spike" in 14 CO 2 caused by atmospheric nuclear weapons testing in the 1940s through 1960s.The two solid curves show model simulations using a firn air transport model either with, or without vertical mixing in the lock-in zone (LIZ).To correctly fit the firn air data, the models used in Buizert et al. (2012) require a LIZ mixing term of the order of 2 × 10 −9 m 2 s −1 ; at the time that study was published, the nature of this mixing remained unclear, and a good fit to the data can be obtained regardless of whether this diffusivity is implemented as molecular-type, eddy-type, or a combination thereof.However, the remnant LIZ mixing is unlikely to be caused by molecular diffusion because of (1) the absence of LIZ gravitational enrichment, and (2) insights from percolation theory that suggest that in a porous medium the gas diffusivity vanishes abruptly at the percolation threshold (i.e., the lock-in depth) as the pore connectivity decreases (Ghanbarian and Hunt, 2014).Here we show that the LIZ mixing can readily be explained by dispersion occurring in the deep firn that is driven by air movement induced through synoptic-scale surface pressure variations.
Steady viscous fluid flow through a disordered porous medium leads to dispersion (Scheidegger, 1954(Scheidegger, , 1961;;Sahimi et al., 1983).On the microscopic level, this process can be understood as a consequence of the different flow paths available to tracer molecules.Figure 2 schematically compares fluid flow through a microscopically ordered (upper) and disordered (lower) porous medium.In the upper panel all pathways are equivalent, and tracer molecules transiting the porous medium each have the same transit time and are not dispersed in passage.In the lower panel the tracer molecules traveling through different sections of the pore space have different transit times, and they are spatially dispersed as they transit the medium.The different pathways not only have unequal lengths.According to the Hagen-Poiseuille equation, the hydraulic conductance of a capillary (for laminar flow) scales as ∝ r 4 , with r the pore radius.Given that a range of pore radii exist in natural firn, the r 4 power law will concentrate flow in the widest pathways, further broadening the distribution of tracer transit times.Figure 1.The 14 CO 2 bomb spike as observed and simulated in firn air at NEEM, northern Greenland (Buizert et al., 2012).Continued mixing within the lock-in zone (LIZ) is needed to correctly simulate the smoothness of the observations.Both simulations use the firn diffusivity values as presented in Buizert et al. (2012); the eddytype LIZ diffusion was either turned on (red curve) or off (blue curve).The case without LIZ mixing was not tuned separately to the observations.On the macroscopic level, dispersion can be described as a diffusive process with a diffusivity D disp with units of m 2 s −1 (Sahimi et al., 1983).Dispersion is much stronger in the longitudinal (along-flow) direction than in the transverse direction.This distinction is important in many applications (for example a point-source contamination in an aquifer subject to groundwater flow), but in firn air applications the source (the atmosphere) is laterally uniform, and only the longitudinal dispersion matters.Fluid flow at a macroscopic mean velocity v results in dispersion of a magnitude where α L is the so-called longitudinal dispersivity, which is a property of the medium with units of m.In the case of firn air dispersion it is important to point out that D disp does not lead to gravitational enrichment, as it originates in macroscopic air movement that does not discriminate between N 2 isotopologues.There are two sources of macroscopic air movement in deep firn: the aforementioned back flux due to pore compaction, and flow driven by synoptic-scale surface pressure variations.Here we will show that the latter is several orders of magnitude larger, and hence dominates the dispersive mixing.The upper few meters of the firn are furthermore subject to wind-driven air flow (e.g.Colbeck, 1989); the dispersion caused by this flow is minimal owing to the low dispersivity of near-surface firn, and overwhelmed by diffusion and convective mixing.
In this work we shall use published gas chromatography experiments by Schwander et al. (1988), which can be used to derive both the diffusivity and the dispersivity of firn samples.Schwander et al. (1988) correctly inferred that molecular diffusion is by far the dominant process in mixing firn air, and that study therefore focused on the first-order process of diffusive transport.We want to emphasize that Jakob Schwander fully realized the potential of these measurements to estimate dispersive gas mixing in firn also, as he wrote: "This dependence of D on the gas-flow velocity can be used to estimate how air flow in the firn caused by atmospheric pressure variations contributes to the gas mixing rate" (Schwander et al., 1988).Since 1988 there has been enormous progress in observing and modeling firn air movement, and presently it is pertinent to also investigate the second-order effect of firn dispersion.Still, the present work should be considered a footnote to the seminal achievement of Schwander et al. (1988), and the many deep and early insights in firn air movement it contains.The reader is strongly encouraged to (re-)read the original work.
2 Air pressure dynamics in the firn column

Mathematical description
Here we present a mathematical description of air pressure dynamics in polar firn, aimed towards understanding firn air movement in deep firn in response to surface pressure variations.
In hydrostatic equilibrium the firn air pressure p(z) increases with depth z following the barometric equation: with p 0 the surface pressure, M the molar mass of air, g the gravitational acceleration, R the gas constant, and T the absolute temperature.Deviations from this hydrostatic balance result in viscous air flow, with a velocity given by Darcy's law: with k the permeability of the firn and µ the (temperatureand pressure-dependent) dynamic viscosity of air.
where k * = k • s op .Equation ( 6) is a nonlinear differential equation that is difficult to solve numerically due to the various terms that scale as O(p 2 ), a complication that arises from the compressibility of air, which makes the fluid density ρ a function of pressure.Equation ( 6) can be simplified substantially by making the following approximation.In the continuity equation Eq. ( 5) the time-varying pressure p(t) is substituted by the mean annual pressure p to get which in turn simplifies Eq. ( 6) to Equation ( 8) is a linear second-order partial differential equation that can be solved easily using finite difference methods.
The approximation in Eq. ( 7) makes the fluid density independent of the pressure; this changes the firn air mass flux from J = ρv = pMv/RT to J = ρv = pMv/RT .By doing so an error is introduced because p(t) can differ from p.This error, however, is rather small given that the atmospheric pressure p never deviates by more than 5 % from the mean annual pressure p.Note that this 5 % error is much smaller than other errors such as the uncertainty in the firn permeability k.Moreover, in the next section it is shown that the pressure response time of the firn is orders of magnitude faster than the synoptic-scale pressure variations that drive barometric pumping in the deep firn, rendering a 5 % error irrelevant.We also neglect the small viscosity changes due to temperature and pressure.

Propagating pressure anomalies into the firn
In the numerical solutions presented here, Eq. ( 8) is solved using a Crank-Nicolson finite difference method, which employs implicit time stepping.At the surface boundary the pressure is set to equal the atmospheric pressure forcing; at the lower boundary the pressure gradient is set to zero.The s op parameterization of Mitchell et al. (2015) is used as well as the k parameterization of Adolph and Albert (2014).
In a first experiment we force the model with an atmospheric pressure p(t) = p + u(t), where u(t) denotes a 1 mbar magnitude unit step (or Heaviside) function at time t = 0. Firn properties of the West Antarctic Ice Sheet (WAIS) Divide site in West Antarctica are used.Equation ( 8) is solved at 0.1 m depth resolution and at time increment t = 1 × 10 −10 year = 3.2 × 10 −3 s.The propagation of the pressure anomaly into the firn is shown in Fig. 3a at four depth levels; the depths z = 67 m and z = 76 m are selected because they are the lock-in depth and deepest firn air sampling depth, respectively.Note that the pressure variations at depth are amplified relative to the surface forcing due to the hydrostatic effect (Eq. 2).Let t1 / 2 (z) be the time required at depth z for the pressure increase to reach half the final amplitude.The depth profile of t1 / 2 (z) is shown in Fig. 4a (note the logarithmic scale).It is clear that pressure fluctuations propagate relatively fast through the firn column; at z = 20 m depth the response time t1 / 2 = 20 s, and at the lock-in depth t1 / 2 = 230 s.This fast response need not be a surprise, given that in free air pressure variations propagate at the speed of sound; in firn, the propagation speed is limited by the finite permeability of the medium.The largest increase in t1 / 2 (z) is seen within the lock-in zone, where firn permeability becomes vanishingly small as s op approaches zero.
The response curves of Fig. 3a can be used via Fourier transform to derive the low-pass filter characteristics of the firn to pressure variations; this is shown in Fig. 3b.For a surface pressure oscillation at any given frequency, the transfer function shows by how much the amplitude of that signal is attenuated in the firn (so a transfer function value of 1 means that the pressure oscillation is transmitted at full amplitude).The model shows that a surface pressure oscillation with a period of 1 h or longer will have a nearly unattenuated response at the LID (yellow curve); pressure variations with a period of 1 min or shorter (such as wind pumping events driven by wind gusts over surface topography) are completely dampened before they reach the LID.These conclusions depend on the assumption that the permeability measured on centimeter-scale samples is also valid on the macroscopic scale of the entire firn column.The ultimate test of these model results would be an in situ pressure gauge buried in the deep firn near the firn-ice transition, but to our knowledge such measurements have only been performed in the upper few meters of snowpacks (e.g.Drake et al., 2016).
The firn transfer function of Fig. 3b can be compared to the spectral density of pressure variability at the site.Figure 3c shows hourly pressure data from an automated weather station (AWS) located at the WAIS Divide site (Lazzara et al., 2012), as well as 6-hourly data from the ERA-Interim reanalysis (Dee et al., 2011); Fig. 3d shows the power spectral densities of both time series.From analyzing the firn transfer functions, it is clear that most of the synoptic-scale pressure variability observed in these records will be expressed at full magnitude all the way down to the firn-ice transition.

Barometric pumping
The pressure variations described above are associated with net macroscopic air movement in the firn column, a phenomenon known as barometric or atmospheric pumping (Nilson et al., 1991).Whenever surface pressure increases, firn air will move downwards as the underlying air parcels are being compressed (the same number of molecules occupy a smaller volume at higher p), resulting in atmospheric air entering the upper firn column.Also, vice versa, when surface pressure decreases, the firn air moves upwards in the column with the upper air parcels being expelled into the free atmo-sphere.One can think of the entire firn column as "breathing" in and out in response to high-and low-pressure systems, respectively.
The 2σ barometrically driven vertical displacement of air parcels ( z) at the WAIS Divide is plotted in blue in Fig. 4b.Near the surface, barometric pumping displaces air vertically of the order of 1 m, which presumably contributes to the establishment of a well-mixed zone (or "convective zone").Note that the curve shows the 2σ vertical displacement; the peak-to-peak displacement is much larger and of the order of 3.5 m near the surface.Figure 4b furthermore shows the mean (absolute) velocity |v| in red.Both curves of course closely resemble each other, given that the vertical velocity is the first derivative of the air displacement.We plotted the actual air velocities in the pores themselves, i.e., the air velocity averaged over the pore cross section.To obtain the velocity per unit of bulk sample cross section, one has to multiply the values by s op .
Another source of macroscopic air movement in deep firn is the gradual closure of the pore space by the densification process, which leads to an upward air flow relative to the firn matrix (Rommelaere et al., 1997).The velocity of this (accumulation-rate-dependent) back flux is of the order of 10 −9 to 10 −8 m s −1 , and is clearly negligible in magnitude compared to the barometrically driven flow.

Experimental dispersivity of polar firn samples
Here we revisit the published firn diffusion experiments by Schwander et al. (1988) to estimate the dispersivity of polar firn.In these experiments, cylindrical firn samples (30 mm diameter, 50 mm length) from Siple station, Antarctica, were placed in a −20 • C chamber inside a gas chromatograph (GC) and flushed through with a pure N 2 carrier gas at a controlled superficial flow velocity ranging from 0.175 to 0.789 mm s −1 .A small amount of CO 2 or O 2 is injected into the N 2 carrier-gas stream before the firn sample, and the eluting peak is measured using a thermal conductivity detector.The effective firn diffusivity (as a function of flow speed) is reconstructed from the peak width in the chromatogram using standard GC theory.More procedural details are given in the original publication (Schwander et al., 1988).
The results at four representative sampling depths are shown in Fig. 5a, in each case at five different flow speeds.The solid lines give the linear least-squares fit to the diffusivity data (second-order terms are small, and neglected here).The (extrapolated) intercept with the y axis (zero flow speed) gives the molecular diffusivity of the sample.The data show that firn diffusivity decreases with depth as the tortuosity of the pore space increases; this result has since been confirmed many times both via direct measurements (Fabre et al., 2000; C. Buizert  Adolph and Albert, 2014) and via inverse methods (Rommelaere et al., 1997;Trudinger et al., 1997;Witrant et al., 2012).
The slope of the fit represents α L , as per Eq. ( 1), with the longitudinal (or along-flow) orientation being vertical in the natural firn setting.Near the surface (z = 2 m) the reconstructed diffusivity is independent of the flow speed, suggesting α L = 0 m.This is due to the homogeneous and open pore geometry of shallow firn, in which the traveled path length and flow speed is identical for all pore clusters that contribute to the flow, and hence tracer molecules traveling through various parts of the pore space do not get dispersed.For deeper firn samples we see an increase in the slope of the fit, and therefore an increase in α L .This increased dispersivity is caused by heterogeneity of the pore geometry, by which the different pore clusters have increasingly disparate path lengths, flow speeds, and cul-de-sacs, which serve to disperse tracer molecules traveling through various parts of the pore space.
The firn dispersivity data are plotted in Fig. 5b as a function of s op (white dots), where we used the average value of the CO 2 and O 2 experiments in Schwander et al. (1988).The data suggest an exponential dependence on s op of the form Propagating the 95 % confidence intervals on the α L values estimated in the slope fitting, the following values for the fitting parameters are suggested: a = 1.26 ± 0.40 m and b = −25.7 ± 2.2.Unfortunately there are no experimental data within the lock-in zone, which makes the exponential extrapolation of Eq. ( 9) uncertain.Using a Monte Carlo scheme, we construct an uncertainty envelope by fitting functions of the form a•exp(b•s op +c•s 2 op ), while forcing the curve to intercept the y axis at an α L (s op = 0) value randomly selected from the interval 0.3 to 3 m, and randomly perturbing the α L data within their 95 % confidence intervals; this envelope is indicated as the blue shading in Fig. 5b.

Dispersive mixing driven by barometric pumping
Multiplying the experimental α L values by the estimated mean absolute velocity profile allows for a theoretical estimate of dispersive mixing in the firn column, as is shown in Fig. 5c for NEEM as the blue curve with the uncertainty envelope.Note that the pore velocities of (Fig. 4b) are first multiplied by s op to convert them into superficial velocities (average velocity in an open tube with the same diameter as the sample) as used by Schwander et al. (1988).The NEEM site is used here because it was the venue for an intensive study of firn air processes by eight different research groups and had a very clear signal of continued mixing within the lock-in zone (Fig. 1).Six different firn air transport models were applied to two separate boreholes at the NEEM site, which together suggest a lock-in zone diffusivity in the range of ∼ 2 × 10 −9 to 7 × 10 −9 m 2 s −1 (Buizert et al., 2012); the theoretical estimate of dispersive mixing derived here corresponds very well with this range.
The dispersivity α L and the pumping velocity |v| have opposite trends with s op .However, because α L scales more strongly with s op , it ends up dominating the behavior of D disp , which generally increases as s op gets smaller -the exception to this pattern is the lower half of the lock-in zone (s op < 0.05), where dispersive mixing decreases again as the air movement approaches zero.Dispersive mixing occurs throughout the firn column; however, within the diffusive zone it is overwhelmed by molecular diffusion, which is about 2-3 orders of magnitude larger.Molecular diffusion effectively ceases at the lock-in depth, and therefore dispersion dominates only in the lock-in depth.

Dispersive mixing in firn air transport modeling
In this section we use the theoretical estimates of dispersive mixing strength in a firn air model to investigate whether it is consistent with the measured trace gas concentrations in air samples extracted from the pore space.We use the NEEM (Buizert et al., 2012) and WAIS Divide (Battle et al., 2011) sites, which are among the most well-characterized firn air sites in the literature.At both locations we use ERA-Interim 6-hourly surface pressure data at the nearest grid point (0.75 • × 0.75 • resolution grid) for the calendar years 2010 and 2011 to calculate the barometrically driven air velocity and dispersive mixing strength.The theoretical D disp values at WAIS Divide are about 15 % higher than those at NEEM (the latter are shown as the blue enveloped curve in Fig. 5c).We assume that dispersive mixing at all times was equal to the average of the 2010-2011 period.
Next, using the CIC firn air model (Buizert et al., 2012), we calibrate the molecular diffusivity profile at both sites using well-established methods (Trudinger et al., 1997;Rommelaere et al., 1997;Battle et al., 2011;Buizert et al., 2012;Witrant et al., 2012;Trudinger et al., 2013), where dispersive mixing strength is set to D disp = γ • D 0 disp , with D 0 disp the theoretical dispersive mixing strength, and γ a scaling factor that is changed in the calibration routine in order to optimize the model fit to the firn air data.The D disp profiles that optimize the fit are shown in Fig. 5c as dashed lines.The fit to the firn air data is optimized at NEEM and WAIS Divide when we set dispersive mixing in the model to around 50 and 58 % of the theoretical estimate, respectively.
The firn air data indicate that WAIS Divide has more dispersive mixing than NEEM, as also predicted by our theoretical calculations.This should thus be considered a robust result.While the theoretical estimates are of the correct order of magnitude, they appear to overestimate the dispersion suggested by observed trace gas concentrations.There may be several causes for this mismatch.First, to fit the same tracer data, different firn air transport models require slightly different diffusivity profiles (Buizert et al., 2012), and some of the mismatch could be an artifact of the firn air model used.Second, at the low s op values of the lock-in zone (where dispersion dominates the transport) our α L parameterization relies on extrapolating the observational estimates (Fig. 5b), and Eq. ( 9) may thus overestimate the true firn dispersivity.Third, the data by Schwander et al. (1988) give the dispersion under steady flow conditions, whereas barometric pumping drives an air flow that is variable in time; the analogy from electronics would be direct current (DC) and alternating current (AC), respectively.Conceivably the AC dispersivity of firn is lower than the DC dispersivity we derived from Jakob Schwander's data, for example because in an alternating flow a portion of the air will retrace its path back to its original position (reducing the dispersion), which does not occur in direct flow.

Gravitational enrichment at Law Dome
A site with strong barometric pumping is the Law Dome site in coastal eastern Antarctica, where firn air has been sampled at the high-accumulation DE08 site (1.2 m year −1 ice equivalent; Etheridge et al., 1996;Trudinger et al., 1997), and the DSSW20K site (0.16 m year −1 ice equivalent; Sturrock et al., 2002;Trudinger et al., 2002).We calculate dispersion at DE08 to be about 65 % stronger than it is at NEEM.
A remarkable property of both Law Dome sites is that LIZ gravitational δ 15 N enrichment is much less than would be expected based on diffusive equilibrium.For DSSW20K and DE08 this is shown in Fig. 6c and e, respectively; note that the shallowest δ 15 N data may be impacted by thermal diffusion due to temperature seasonality (Severinghaus et al., 2001).Here we will focus on the DSSW20K site, for which more data are available.At DSSW20K gravitational enrichment in δ 15 N appears to stop ∼ 7 m above the actual lockin depth.Here we explore the possibility that this anomaly is due to dispersive mixing in the deep firn.First, we calibrate the CIC firn model to the DSSW20K site using firn air data of CH 4 , CO 2 , SF 6 , 14 CO 2 , CFC-11, CFC-12, and CFC-113 (Sturrock et al., 2002) using established methods (Buizert et al., 2012).Following our findings at the NEEM and WAIS Divide sites, we set the deep firn dispersion to equal 55 % of our theoretical estimate based on Law Dome surface pressure time series from the ERA-Interim reanalysis.Next, we explore four instructive modeling scenarios that are color-coded in Fig. 6.
In the first scenario (blue curves), we have eliminated both near-surface convection and deep-firn dispersion to show the δ 15 N gravitational signal in the absence of macroscopic mixing (Fig. 6c).It is clear that the δ 15 N data in the lock-in zone are depleted about 40 per meg relative to gravitational equilibrium.In the second scenario (red curves) we add the convective and dispersive mixing, and find an improved fit to all tracers except δ 15 N.This second scenario is equivalent to best current modeling practices, and comparable to methods used above at the NEEM and WAIS Divide sites.We must conclude that our best practices cannot account for the anomalous deep-firn δ 15 N signal seen at DSSW20K (Fig. 6c) and DE08 (not shown).
In the third scenario (yellow curves), we attempt to fit the δ 15 N data simply by increasing the magnitude of D disp .We have to increase dispersion 150-fold to simulate δ 15 N correctly; however doing so seriously compromises the model fit to all other tracers (Fig. 6a-b), showing this approach to be invalid.
In all modeling scenarios so far we have assumed that molecular diffusion and dispersion are both one-dimensional processes that vary smoothly with depth and that occur inde- pendently without interactions between them.In reality, the pore space is a three-dimensional network that is strongly impacted by density layering (Martinerie et al., 1992;Hörhold et al., 2011;Mitchell et al., 2015).Diffusion and dispersion have generally opposite relationships to density.Thus, near the lock-in zone of a firn with strong density layering, we should expect to see alternating bands of high diffusivity (associated with low-density strata) and bands with high dispersion (associated with high-density strata).The high-density strata have a larger fraction of closed bubbles, and therefore relatively few vertically connective pathways; these pathways will channel the barometrically induced flow, thereby becoming focal points of dispersive mixing.This situation is depicted schematically in Fig. 7a.In this conceptual model, dispersion dominates transport in the vertical direction, which leads to a strongly reduced gravitational enrichment.However, in the transverse directions molecular diffusion still dominates, particularly in the low-density strata.
In a layered firn it is thus conceivable that a mixed zone exists between the molecular-diffusion-dominated diffusive zone and the dispersion-dominated lock-in zone, in which gravitational enrichment is very weak, yet molecular diffusion is still active.
Next we attempt to capture the dynamics of such a layered firn in our one-dimensional model.In our fourth modeling scenario (green curves), we use an idealized layered firn model, with alternating annual bands of diffusive and dispersive mixing; details are shown in Fig. 7.The dispersion has been set to perfectly compensate the reduced molecular diffusion, and therefore the total mixing is comparable to that of scenario 2. We let the compensation be perfect for CO 2 , which means it is imperfect for other trace gases.We find that using such a layered approach can simulate both the regular tracers and δ 15 N correctly (Fig. 6a-c).These modeling results confirm that δ 15 N enrichment can cease before the lock-in depth is reached.The common practice of defining the lock-in as the depth where δ 15 N enrichment stops (e.g., Sowers et al., 1992) may thus be invalid at sites that have both strong layering and large barometric variations.Instead, the lock-in depth should be identified with transient tracers such as CO 2 and CH 4 , using their sharp inflection point.However, our fourth modeling scenario is still flawed as it attempts to represent what is fundamentally a three-dimensional process into a one-dimensional model.Furthermore our choice of perfectly compensating variations in dispersive and diffu- The Cryosphere, 10, 2099-2111, 2016 www.the-cryosphere.net/10/2099/2016/sive mixing is of course questionable -it is included here as an illustrative example only.To gain a meaningful representation of these processes it may be necessary to move to firn air transport models of two or more dimensions.

Implications for firn air and ice core studies
There are probably three factors that contribute to the magnitude of dispersive mixing at any given site.
1. Magnitude of barometric variability.Figure 8 shows maps for Greenland and Antarctica of the time-averaged root mean square of the surface pressure change dp/dt expressed in mbar day −1 ; this is a good proxy for the barometric pumping power available at a given site, as the air flow velocity scales with the rate of pressure change.The maps suggest that barometric pumping is strongest near the ice sheet margins, and weakest in the interior.The Dome A ice core drilling site at Kunlun station has the smallest barometric pumping of all sites, whereas coastal cores such as Law Dome and James Ross Island (JRI) have strong barometric pumping.
2. Firn column thickness.The amount of barometrically driven air flow at any given depth depends on the total size of the air reservoir below that depth (which is why air flow decreases with depth in Fig. 4).All other things being equal, a thicker firn column will have higher air flow velocities and hence more dispersive mixing than a thin column.

3.
Layering.As argued in Sect.3.4, firn density layering can enhance dispersive mixing by increasing medium heterogeneity, specifically through the formation of high-density, high-dispersivity layers.Melt layers and ice lenses may similarly act as focal points for dispersive mixing.
The NEEM and WAIS Divide sites have comparable firn thickness and density layering, and therefore the stronger barometric variability at the WAIS Divide site results in stronger dispersive mixing at that site (Fig. 5).The Law Dome DE08 and DSSW20K sites experience the same barometric variations, yet DE08 has a thicker firn column and DSSW20K has more (annually spaced) high-density layers.It is therefore not clear a priori which of these two Law Dome sites has stronger dispersive mixing, and unfortunately the available firn air data are of insufficient resolution to establish this unambiguously.
Dispersive mixing influences the ice core record in several ways, the most important of which is via the broadening of the gas age distribution.A comparative firn model study at the NEEM site showed that the low-diffusion lock-in zone environment contributes more to the broadening of the final age distribution than the diffusive zone does (Buizert et al., 2012).The weak barometric variability in interior Antarctica (Fig. 8) may be part of the reason why gas records from central East Antarctic ice cores such as EPICA Dome C have surprisingly little smoothing, as evidenced by the fact that abrupt methane transitions during e.g., the Younger Dryas period and Dansgaard-Oeschger cycle are well preserved (Loulergue et al., 2008).
Dispersive mixing potentially has implications for the use of δ 15 N as a proxy for past firn column thickness, which is an additional constraint on past age (Sowers et al., 1992;Schwander et al., 1997;Goujon et al., 2003;Parrenin et al., 2012;Capron et al., 2013;Buizert et al., 2015).The cited studies all assume δ 15 N reflects the thickness of the diffusive zone, which is then used to estimate the lock-in depth by adding an estimated convective zone thickness.As we showed in Sect.3.4, under circumstances of strong layering and intense barometric pumping, δ 15 N may underestimate the lock-in depth where age is fixed -this effect is likely to be important at coastal coring sites such as Law Dome, James Ross Island, Berkner Island, and Roosevelt Island where barometric variability is strong, as well as at sites influenced by melt layers.It has been suggested that glacial firn has more pronounced layering than present-day firn (Bendel et al., 2013).However, climate models participating in PMIP2 (Paleoclimate Modeling Intercomparison Project Phase 2) disagree on the sign and magnitude of the change in cyclonic activity around Antarctica between the Preindustrial and Last Glacial Maximum (Rojas et al., 2009).It is therefore conceivable, but highly uncertain, that dispersive mixing was stronger during glacial times.Several studies have noted a δ 15 N model-data mismatch in central Antarctica during glacial climate conditions, with densification models simulating a thickening of the firn column (relative to present) and δ 15 N data suggesting a thinning (Landais et al., 2006;Capron et al., 2013).We speculate that enhanced C. Buizert and J. P. Severinghaus: Dispersion in firn by pressure variability glacial dispersive mixing could contribute to the observed low gravitational enrichment in δ 15 N. Overall, more multitracer, high-resolution firn air studies at sites with strong barometric variability will be needed to better understand the influence of barometric pumping on gravitational enrichment and the δ 15 N proxy.

86 Kr excess as a potential proxy of past synoptic activity
The degree of isotopic gravitational enrichment of any given gas species in the firn depends on the relative strength of molecular diffusion, which acts to drive isotopic enrichment towards gravitational equilibrium, and macroscopic transport processes (convection, advection and dispersion), which act to erase the enrichment.Slow-diffusing gases such as krypton (Kr) and xenon (Xe) will therefore always be less isotopically enriched than fast-diffusing gases such as N 2 and argon (Ar).This effect was first observed by Kawamura et al. (2013), who studied deep air convection (a form of macroscopic transport) at the Antarctic Megadunes site (Severinghaus et al., 2010).
Here we define 86 Kr excess as 86 Kr excess = δ 86 Kr / 82 Kr| tc − δ 40 Ar / 36 Ar| tc , where the subscript "tc" ("thermally corrected") denotes that the isotopic ratios have been corrected for thermal fractionation either by using a thermal model or by paired δ 15 N-δ 40 Ar data (Grachev and Severinghaus, 2003a, b).Due to the different molecular diffusivities of N 2 and Kr, they are in a different state of gravitational disequilibrium, which makes 86 Kr excess a measure of the aggregate strength of non-diffusive transport processes in the firn column.There are several reasons for using 86 Kr and 40 Ar rather than e.g., 136 Xe and 15 N. First, Kr is over 10 times as abundant in the atmosphere as Xe and therefore more easily measured.Second, δ 40 Ar is more precise in terms of gravitational signal-to-noise ratio, and less affected by thermal diffusion than δ 15 N. Third, δ 40 Ar and δ 86 Kr can both be measured on air extracted from the same physical ice sample after removal of O 2 , N 2 , and other reactive gases via gettering (Severinghaus et al., 2003).Because δ 86 Kr is always less gravitationally enriched than δ 40 Ar, 86 Kr excess is always negative.
Figure 6d shows 86 Kr excess for all four modeling scenarios discussed in Sect.3.4.In scenario 1 (blue curves) the only macroscopic transport mechanism is advection, which results in a 86 Kr excess of −5 per meg.Adding convection and modest dispersion (scenario 2, red curve) increases 86 Kr excess in magnitude to −10.2 per meg, reflecting the increased degree of gravitational disequilibrium in the firn column.Both scenarios with strong dispersion (scenarios 3 and 4, yellow and green curves, respectively) show a further increase in 86 Kr excess magnitude to 16.5-18 per meg.
Measurements on the WAIS Divide ice core (WAIS-Divide Project Members, 2013; WAIS Divide Project Mem-bers, 2015) show 86 Kr excess values of around −35 per meg during the late Holocene (analytical precision is better than 20 per meg).However, older sections of the core show 86 Kr excess values as low as −90 per meg (Orsi et al., 2014), suggesting periods of greatly enhanced gravitational disequilibrium in the firn column (data not shown).One tantalizing interpretation could be that these very negative 86 Kr excess values represent periods of enhanced synoptic activity in (West) Antarctica, driving strong dispersive mixing.
We propose here that 86 Kr excess may act as a proxy for past synoptic-scale pressure variability (or storminess/cyclone density).Before this interpretation can be accepted, however, additional work is needed.The large spatial variability in barometric variability over Antarctica (Fig. 8) provides a valuable opportunity to verify whether ice core 86 Kr excess indeed scales with local barometric variability.Additional work is needed to reliably correct 86 Kr excess for the influence of advection and convection in the firn column et al., 2013), preferably through detailed studies of 86 Kr excess in modern-day firn.If corroborated, 86 Kr excess could hold important clues about past changes to the large-scale atmospheric circulation, particularly when combined with reconstructions of past wind conditions from surface roughness (Barnes et al., 2006;Bay et al., 2010).In Greenland, changes in synoptic activity are presumably linked to the buildup of the Laurentide ice sheet and its orographic influence on the storm tracks (Kageyama et al., 1999).In Antarctica, synoptic activity may be linked to e.g., meridional movement of the eddy-driven jet (Ceppi et al., 2013), atmospheric teleconnections to the tropical Pacific (Ding et al., 2012), or the position and depth of the Amundsen Sea low (Turner et al., 2013).

Summary and conclusions
In this work we show that surface pressure variability on synoptic timescales drives macroscopic air movement in the deep firn, which in turn leads to dispersion of trace gases in the firn open porosity.The work resolves an outstanding question regarding the nature of lock-in zone mixing deduced from detailed firn air experiments at the north Greenland NEEM site (Buizert et al., 2012).
We present a mathematical description of the propagation of pressure anomalies in polar firn.We find that pressure variations on the timescale of order 1 h or slower are propagated to the firn-ice transition at full amplitude; variations on shorter timescales are attenuated.Net barometrically driven air movement is on the centimeter scale in the deep firn, and on the meter scale in the upper firn; mean velocities are of the order of 10 −6 m s −1 .The precise values of the air displacement and velocity depend primarily on the firn thickness and the barometric variability at the site.
We use published firn sample gas chromatography experiments to estimate the dispersivity α L of firn, which is the proportionality constant in the relationship between superfi-cial gas velocity and the apparent dispersion strength (Eq.1).We find that α L scales exponentially with the open porosity s op , with α L ≈ 0.1 m at the lock-in depth.Combining simulated air velocities and firn dispersivity, we calculate dispersive mixing in the deep firn to be of the order of 10 −9 m 2 s −1 , with precise values again depending on firn thickness and the barometric variability at the site.
We apply these theoretical estimates of dispersion in a firn air transport model, and find that they overestimate the amount of lock-in zone dispersivity needed to optimize the fit to firn air trace gas measurements; this mismatch may be due to the fact that our firn dispersivity parameterization is based on steady-flow conditions, whereas barometric pumping induces a time-variable flow.We suggest that strong dispersive mixing at Law Dome, Antarctica, in combination with firn layering, may halt gravitational enrichment in δ 15 N before the lock-in zone is reached.
The dispersive mixing discussed here increases scientific understanding of firn air transport, and has direct implications for the modeling thereof.The ice core record is impacted by dispersive mixing primarily through the widening of the gas age distribution.We propose that 86 Kr excess may be an ice core proxy for past synoptic activity, which is linked to the large-scale atmospheric circulation.
6 Data availability NEEM and WAIS Divide firn air data are available with the original publications (Buizert et al., 2012;Battle et al., 2011).Siple station firn sample gas chromatography data are listed in Table 1 of Schwander et al. (1988).ERA-Interim reanalysis products are publicly available through the European Centre for Medium-Range Weather Forecasts (http://www.ecmwf.int).Automated weather station (AWS) data can be requested from the University of Wisconsin-Madison Automatic Weather Station program (http://amrc.ssec.wisc.edu/).Model code and output are available upon request from the first author (buizertc@oregonstate.edu).

Figure 2 .
Figure 2. Schematic depiction of dispersion.Fluid movement through a microscopically homogeneous, low-tortuosity porous medium (upper panel) is non-dispersive; fluid movement through a microscopically disordered porous medium (lower panel) disperses tracer molecules.Macroscopic fluid flow is from left to right.

Figure 3 .
Figure 3. Propagation of pressure variations into the firn at WAIS Divide.(a) Pressure response at indicated depths to a unit-step surface pressure increase of 1 mbar at time t = 0. (b) Low-pass filter characteristics of the firn column to surface pressure variability.(c) Surface pressure history of year 2010 at the WAIS Divide site from AWS measurements (orange) and from the ERA-Interim reanalysis (green).(d) power spectral density (PSD) of the time series shown in (c).

Figure 4 .
Figure 4. Firn air pressure response at the WAIS Divide site.(a) Response time t1 / 2 (z) to a unit-step surface forcing.(b) The 2σ air displacement z as a function of depth (blue), and the mean absolute air velocity due to barometric pumping (red), both forced by the 2010 surface pressure time series (Fig. 3c).

Figure 5 .
Figure 5. Firn dispersivity.(a) Effective diffusivity of CO 2 as a function of flow speed as measured by Schwander et al. (1988) on firn samples from Siple station.(b) Experimental firn longitudinal (i.e., along-flow) dispersivity α L as a function of sample open porosity (white dots) with best fit (solid blue line); the shaded area indicates the estimated 95 % confidence interval of the fit.(c) Dispersive mixing at NEEM, using the dispersivity envelope (95 % confidence) from (b), multiplied by the mean absolute firn air velocity induced by barometric pumping at that site.The dashed black (red) line indicates the optimal level of dispersive mixing needed in the firn air transport model to fit the NEEM (WAIS Divide) firn air observations(Buizert et al., 2012;Battle et al., 2011).

Figure 6 .Figure 7 .
Figure 6.Model simulations and firn air data for the Law Dome DSSW20K (a-d) and DE08 (e) sites, with model scenarios color-coded.The first scenario (CTRL, blue) is a control run with no convective mixing or dispersion to show gravitational equilibrium.The second scenario (+CZ + D disp , red) includes convective mixing (following Kawamura et al. (2006) with H = 2 m and D eddy,0 = 1 × 10 −5 m 2 s −1 ) and dispersion to optimize the fit to all tracers except δ 15 N.The third scenario (150 × D disp , yellow) uses strongly enhanced dispersive mixing to optimize the fit to δ 15 N data.The fourth scenario (layered D disp , green) uses layered dispersion and diffusivity to optimize the fit to all tracer data.Model details and atmospheric trace gas forcings are as inBuizert et al. (2012).Spatial model resolution is z = 0.1 m for scenarios 1-3, and z = 0.01 m for scenario 4. (a) DSSW20K methane mixing ratio; (b) DSSW20K 14 CO 2 mixing ratio; (c) DSSW20K δ 15 N-N 2 ; (d) DSSW20K 86 Kr excess; (e) DE08 δ 15 N-N 2 .

Figure 8 .
Figure 8. Synoptic-scale barometric pumping strength for Greenland and Antarctica using 6-hourly ERA-Interim reanalysis surface pressure values for the period 1 January 2010 through 31 December 2011.The color scale gives the root mean square of the pressure rate of change dp/dt in units of mbar day −1 .Key ice core drilling locations are indicated.