The influence of layering and barometric pumping on firn air transport in a 2-D model

Ancient air trapped in ice core bubbles has been paramount to developing our understanding of past climate and atmospheric composition. Before air bubbles become isolated in ice, the atmospheric signal is altered in the firn column by transport processes such as advection and diffusion. However, the influence of low-permeability layers and barometric pumping (driven by surface pressure variability) on firn air transport is not well understood and is not readily captured in conventional one-dimensional (1-D) firn air models. Here we present a two-dimensional (2-D) trace gas advection–diffusion–dispersion model that accounts for discontinuous horizontal layers of reduced permeability. We find that layering or barometric pumping individually yields too small a reduction in gravitational settling to match observations. In contrast, when both effects are active, the model’s gravitational fractionation is suppressed as observed. Layering focuses airflows in certain regions in the 2-D model, which acts to amplify the dispersive mixing resulting from barometric pumping. Hence, the representation of both factors is needed to obtain a realistic emergence of the lock-in zone. In contrast to expectations, we find that the addition of barometric pumping in the layered 2-D model does not substantially change the differential kinetic fractionation of fastand slow-diffusing trace gases. Like 1-D models, the 2D model substantially underestimates the amount of differential kinetic fractionation seen in actual observations, suggesting that further subgrid-scale processes may be missing in the current generation of firn air transport models. However, we find robust scaling relationships between kinetic isotope fractionation of different noble gas isotope and elemental ratios. These relationships may be used to correct for kinetic fractionation in future high-precision ice core studies and can amount to a bias of up to 0.45 C in noble-gas-based mean ocean temperature reconstructions at WAIS Divide, Antarctica.


Introduction
In the upper 50-130 m of consolidated snow above an ice sheet, known as the firn layer, atmospheric gases gradually become entrapped in occluded pores and are eventually preserved as bubbles in the ice below.Antarctic ice core records containing these trapped gases have been critical in informing our understanding of the interplay of past climate and atmospheric trace gas variability over the past 800 000 years (Petit et al., 1999;Lüthi et al., 2008).As atmospheric gases migrate through the firn, they are modified in elemental composition and isotopic signature by several competing physical processes (Schwander et al., 1988(Schwander et al., , 1993;;Trudinger et al., 1997;Buizert et al., 2012;Kawamura et al., 2013;Mitchell et al., 2015).Therefore, appropriate corrections must be applied to firn and ice core records to accurately reconstruct atmospheric trace gas histories.
Firn is a layered medium, in which the denser layers can impede vertical diffusion and transport (Hörhold et al., 2012;Mitchell et al., 2015;Orsi et al., 2015) (Fig. 1).The significance of these layers for firn gas transport remains unclear and motivates this work.Readers who are familiar with the structure of firn and its air transport processes may wish to skip ahead to the last paragraph of this section.To build some Published by Copernicus Publications on behalf of the European Geosciences Union.intuition about firn transport processes, a commonly used analytical model of firn air transport is provided in Appendix A.

Box 1| Porous media terminology
Porosity: the fraction of (firn) volume filled by gas Permeability: the degree to which a porous medium permits viscous flow to pass through Fickian diffusion: molecular diffusion that is proportional to the concentration gradient as described by Fick's first law Tortuosity: measure of the twistedness of pathways through a porous medium We distinguish four main processes affecting the composition of air in firn: diffusion, advection, dispersion, and convective mixing.Molecular diffusion, driven by concentration gradients in the firn, is the primary mode of horizontal and vertical transport.Molecular diffusion also enables gravitational fractionation, or "settling", of trace gases in proportion to their masses (Schwander, 1989;Sowers et al., 1989;Schwander et al., 1993;Trudinger et al., 1997).Gravitational settling leads to an enrichment of heavy isotopes with depth that is described in equilibrium by the barometric equation (Schwander, 1989;Sowers et al., 1989;Craig et al., 1988): where δ ≡ r sample r standard − 1 ≡ q − 1, with r being the isotope ratio (unitless), z the depth (m), T the absolute temperature (K), m the isotope mass difference (kg mol −1 ), g the gravitational acceleration (m s −2 ), and R the fundamental gas constant (J mol −1 K −1 ).
Gradual accumulation of snow and air bubble trapping leads to a slow, downward advection of the enclosed air.The net air advection velocity is slower than the snow accumulation rate (yet still downward in an Eulerian framework) because compression of the porous firn medium produces a flow of air upward relative to the firn matrix (Rommelaere et al., 1997).
Buoyancy-driven convection and brief pressure anomalies associated with wind blowing over an irregular topography cause strong mixing between the near-surface firn and the unfractionated atmosphere, smoothing concentration gradients (Colbeck, 1989;Severinghaus et al., 2010;Kawamura et al., 2013).This mixing causes substantial deviations from the gravitational settling equilibrium (i.e., the solution to Eq. 1) and leads to varying degrees of kinetic isotope fractionation.This is because faster-diffusing isotopes more readily return to thermal-gravitational equilibrium by diffusion (Buizert et al., 2012;Kawamura et al., 2013).
Lastly, surface barometric pressure variability on longer timescales (> 1 h) drives air movement down to the firn-ice transition.Building on work by Schwander et al. (1988), Buizert and Severinghaus (2016) suggest that surface pressure variability may produce significant pressure gradients in the firn that induce airflow.Porous firn has a high tortuosity, i.e., two points are typically connected by strongly curved paths, and the deep firn also contains many cul-desacs (Buizert and Severinghaus, 2016, their Fig. 2).Airflow through such a medium produces mass-independent dispersive mixing.Dispersion in this context is an emergent macroscopic phenomenon that describes microscopic velocity de- viations from Darcy's law of bulk fluid flow in different pores.This process may be accounted for by adding a dispersive mixing term to the advection-diffusion equation traditionally used for trace gas transport in firn (e.g., Buizert and Severinghaus, 2016).Together, these four processes yield a firn column that is split into a surface mixed zone (SMZ, which has historically been labeled the convective zone), a diffusive zone (DZ), and a lock-in zone (LIZ) (Fig. 2).The close-off depth (COD) occurs where the air content becomes fixed and pressure in open pores increases above hydrostatic pressure.We prefer the term "surface mixed zone" over the more commonly used terminology "convective zone" to acknowledge the dual nature of mixing driven by convection and high-frequency pressure variability in this region.Large seasonal temperature gradients in the SMZ can lead to isotopic fractionation, which is only partially attenuated (Severinghaus et al., 2001(Severinghaus et al., , 2010;;Kawamura et al., 2013).Molecular diffusion dominates in the DZ but effective firn diffusivity decreases with depth due to the increasing influence of tortuosity hindering diffusion.Throughout the DZ, gravitational settling leads to an enrichment of all isotopes heavier than air in proportion to their mass difference.The top of the LIZ, the somewhat poorly defined lock-in depth (LID) horizon, is commonly deduced from a rather sudden change in the slope of the δ 15 N, CO 2 , or CH 4 profiles.Gravitational enrichment of isotopes ceases in the LIZ and isotope ratios remain constant with depth.We term any such deviation from gravitational equilibrium "disequilibrium" (without implying that such a situation is not in steady state).
The physical mechanism responsible for the cessation of gravitational enrichment in the deep firn is still not fully understood.Since chlorofluorocarbons (CFCs) and other anthropogenic tracers have been detected in firn air measurements from the LIZ well below the depth expected from pure advection, it is clear that some amount of vertical transport by molecular diffusion or dispersion continues in the LIZ (Severinghaus et al., 2010;Buizert et al., 2012;Buizert and Severinghaus, 2016).However, no further gravitational settling of isotopes occurs in the LIZ as indicated by constant δ 15 N values.Furthermore, the vertical transport in the LIZ appears to be at least to some degree dependent on mass and diffusivity since the faster-diffusing CH 4 advances further in the LIZ than the slower-diffusing gases CFC-113 or CO 2 (Buizert et al., 2012).Therefore, transport in the LIZ cannot be explained by either mass-indiscriminate dispersive mixing or molecular diffusion alone.Most current 1-D firn air models use a greatly reduced molecular diffusivity in the LIZ and simultaneously introduce a mass-independent mixing term tuned to match measured trace gas profiles (Buizert et al., 2012).
Here, we explore the possibility that non-fractionating trace gas mixing in deep firn may be explained by the combination of barometric pumping and discontinuous horizontal layers that have nominal diffusivity.High-density layers are empirically linked to low vertical permeability, increasing the firn's tortuosity and forcing extensive horizontal transport.The influence of layering on firn gas transport is mostly untested in numerical models so far since previous firn air models were generally limited to one dimension.In particular, we will test two mechanisms by which density layering could influence isotope ratios in firn air: (i) layering may reduce gravitational settling of isotopes because vertical settling of isotopes is absent during horizontal transport along layers; (ii) layering may modulate the mass-independent dispersive mixing effect of barometric pumping.Our analyses will focus on two Antarctic high-accumulation sites, WAIS Divide and Law Dome DSSW20K (Trudinger et al., 2002;Battle et al., 2011).

Governing equation and firn properties
We model 2-D trace gas transport in firn by numerically solving the advection-diffusion-dispersion equation, known from hydrology (Freeze and Cherry, 1979), adapted to firn (following Schwander et al., 1993;Rommelaere et al., 1997;Trudinger et al., 1997;Severinghaus et al., 2010;Buizert et al., 2012;Kawamura et al., 2013;Buizert and Severinghaus, 2016): with q ≡ δ + 1 the ratio of any isotope to 28 N 2 relative to a standard material, s ≡ s op exp m g R T z the pressurecorrected open porosity (m 3 m −3 ), T temperature (K), the thermal diffusion sensitivity (K −1 ), and u the advection velocity field (m s −1 ).D m and D d are the 2-D molecular diffusion and dispersion tensors (m 2 s −1 ).∇ q is the concentration gradient and ∇• denotes the 2-D divergence operator.From left to right, the terms of Eq. ( 2) represent the rate of change in concentration or isotope ratio, Fickian diffusion, gravitational fractionation, thermal fractionation, dispersive mixing, and advection.Since Eq. ( 2) is only valid for the binary diffusion of a trace gas into a major gas, ratios of any two isotopes of masses x and y are obtained by separately simulating the transport of each isotope into the major gas 28 N 2 and using the relationship q x/y = q x/28 q y/28 (3) to calculate the isotope ratios of interest (Severinghaus et al., 2010).Isotope ratios are assumed to be constant at the surface (Dirichlet boundary) and reconstructions of atmospheric  to force runs of these anthropogenic tracers (see Supplement).The model time step is 3.5 days; smaller time steps have little impact on the model results and make the model impractical to run due to computational costs.The bottom boundary only allows for the advective flux to leave the domain (Neumann boundary).Diffusion and dispersive mixing cease below the COD.A periodic boundary condition is used in the horizontal direction.The horizontal extent of the model is varied between sites with differing snow accumulation rates to maintain a constant ratio of annual layer thickness to the model's spatial extent, which affects barometric pumping velocity.Firn density (Fig. 3a) is prescribed from a fit to the measured density profile at each site.Following Severinghaus et al. (2010) and Kawamura et al. (2013), empirical relationships are used to derive open and closed porosities from the density profile (Fig. 3b).The pressurecorrected open porosity s is assumed to be time independent.

Advection velocity and barometric pumping
The 2-D velocity field u is a result of a combination of (i) air migration with the firn (w firn ), (ii) return flow of air from the firn to the atmosphere due to the gradual compression of pores (u r ), and (iii) airflow resulting from barometric pumping (u b ) (Figs. 3c, 4).w firn is constrained by assuming a timeconstant snow and ice mass flux at all depths.u r is calculated based on the effective export flux of air in open and closed pores at the COD, imposing a constant mean vertical air flux throughout the firn column (Fig. 3d) (Rommelaere et al., 1997;Severinghaus and Battle, 2006).u b is the airflow needed to re-establish hydrostatic balance in the firn in response to any surface pressure anomaly.Surface pressure variability is represented by (pseudo-) red noise, mimicking observed pressure variability at both sites.The near-coast location Law Dome is more strongly affected by storm activity than WAIS Divide with pressure variability 11.2 hPa day −1 compared to 7 hP day −1 at WAIS Divide.u r and u b follow Darcy's law of flow through porous media (Darcy, 1856): with ∇ P the pressure gradient, κ the permeability of firn, and µ the viscosity of air (Fig. 4).Further details on the derivation of these velocities are provided in the Supplement.

Firn layering
Idealized firn layering is implemented by forcing the vertical velocity components u r and u b , as well as all vertical diffusive fluxes between the grid boxes on either side of a layer to be zero.Only the advection of air with the firn (w firn ) remains active at these grid box boundaries.Layering limits vertical gas migration and yields almost exclusively horizontal transport between layers.We assign layers an infinitesimal thickness because of the computationally limited spatial resolution of the model.Layers are repeatedly introduced at a specific depth and migrate down with the velocity of the firn.The vertical distance between layers is set to correspond to the snow accumulation of 1 year and the horizontal extent of layers increases linearly with depth until they cover the entire domain at the COD.The mean layer opening size is held proportional to the annual layer thickness to make the vertical advection velocities independent of the arbitrary horizontal extent of the model.To obtain more realistic flow fields, the permeability of layers increases gradually towards both ends of a layer.

Dispersive mixing
The dispersion tensor D d is made up of two components: (i) non-fractionating mixing of air in the SMZ and (ii) dispersive mixing caused by barometric pumping in the tortuous firn medium.First, the SMZ is commonly represented by mass-independent ("eddy") diffusion acting in the vertical and horizontal directions.The corresponding diffusivity profile D SMZ (z) is prescribed as an exponential decay away from the snow-atmosphere interface (Kawamura et al., 2013).Its maximum surface value and the decay constant are chosen to match observed δ 15 N values in the deep firn.After reaching a specified maximum depth of 8 m at WAIS Divide and 14 m at Law Dome DSSW20K, D SMZ tapers linearly to zero over the next 2 m.Second, airflow through any dispersive medium leads to mixing in the directions longitudinal and transverse to the flow.Because barometric pumping velocities are orders of magnitude faster than the return flow, dispersive mixing is dominated by barometric pumping.The degree of dispersive mixing in firn presumably depends on the direction of flow and differs between the longitudinal-to-flow and transverseto-flow direction.However, the treatment of anisotropic media is complex and only one parameterization for vertical, longitudinal-to-flow dispersion in firn is currently available (Buizert and Severinghaus, 2016).Therefore, we assume that the dispersivity α of firn is isotropic and linearly dependent on the magnitude of the flow velocity vector (v ≡ |u b + u r |).In this case, the 2-D dispersion tensor becomes with I the second-order identity matrix (Freeze and Cherry, 1979;Buizert and Severinghaus, 2016).The dispersion flux term in Eq. ( 2) simplifies to The dispersivity parameterization of Buizert and Severinghaus ( 2016) is based on direct measurements of cylindrical firn samples from Siple Station, Antarctica, performed by Schwander et al. (1988).The parameterization relates dispersivity to open porosity s op as A factor of s was added to the original parameterization by Buizert and Severinghaus (2016) because α relates dispersive mixing to flow velocities per unit pore cross section (w pores ).Schwander et al. (1988), however, originally measured the considerably slower bulk airflow per unit firn cross section (i.e., w bulk = w pores s ).Since dispersivity is a scale-dependent property, it is important to use parameterizations that are compatible with the resolution of the numerical model.The sample size of Schwander et al. (1988) (30 mm diameter and 50 mm length) approximately matches the resolution of our numerical model (30 × 40 mm) and thus should adequately approximate subgrid-scale (i.e., pore-scale) mixing processes that currently cannot be resolved.Spatial inhomogeneity of subgrid-scale firn dispersivity that was not captured by the sampling of Schwander et al. (1988) cannot be accounted for in the model.Dispersion on larger scales such as the interaction of flow and layers is explicitly represented in the model by the interplay of advection and diffusion.Thus, dispersive mixing is fully constrained in the model and based on empirical parameterizations that are not subject to tuning.

Molecular diffusion
The (effective) molecular diffusivity profile is established by simultaneously fitting the simulated CO 2 and CH 4 profile to real firn measurements at both sites.Effective vertical diffusivity decreases with depth to represent the subgridscale effect of decreasing pore connectivity and increasing firn tortuosity, which is not fully represented by the explicit macroscopic layers in our model (Fig. 5).Diffusivities for other trace gases are calculated by scaling the tuned CO 2 diffusivity by the free air diffusivity of each gas relative to CO 2 (Trudinger et al., 1997).The diffusivity tuning presents an underconstrained problem because horizontal and vertical molecular diffusivities are essentially free parameters.It is qualitatively evident from firn air sampling that horizontal connectivity or diffusivity is much higher than vertical diffusivity in deep firn, but no satisfactory quantification of this anisotropy is available.As a best guess estimate, we set horizontal molecular diffusivities equal to 10× the vertical diffusivity at the same depth.There are many degrees of freedom in tuning molecular diffusivities and our diffusivity parameterization is therefore not unique.However, sensitivity tests with equal horizontal and vertical diffusivity in the model (compensated by shorter horizontal layers) yield comparable results.Considerable temperature gradients can exist in presentday firn because of recent global atmospheric warming and these gradients can lead to increased isotope thermal fractionation, in particular of δ 15 N.The sensitivity of isotopes to diffuse in response to temperature gradients is captured by the thermal diffusion sensitivity .The temperature dependence of is approximated as a function of the effective average temperature T in Kelvin: or is assumed to be temperature independent if the coefficients a and b are unknown for a specific isotope ratio (Severinghaus et al., 2001).Coefficients a and b were determined experimentally for different isotope ratios by Grachev and Severinghaus (2003a, b) and Kawamura et al. (2013).
A table of all model parameters and further details of the numerical realization of 2-D gas transport is provided in the Supplement.

CO 2 and CH 4
A comparison of simulated and observed CO 2 and CH 4 profiles shows good agreement at WAIS (Fig. 6).In line with observations, both CO 2 and CH 4 concentrations decrease slowly with depth until ∼ 68 m below the surface due to the gradually increasing gas age and the anthropogenic rise in atmospheric CO 2 and CH 4 concentrations.The more rapid decrease of CO 2 and CH 4 below ∼ 68 m is explained by a much slower vertical penetration of air and a faster increase of the gas age with depth in the LIZ.
In the following discussion we will examine and compare results from four different versions of the 2-D model: with or without impermeable layers and with activated or deactivated barometric pumping.In versions without layering, our 2-D model loses all horizontal heterogeneity and will thus be referred to as a "1-D model" from here on.Since the explicitly implemented tortuosity from layering in the 2-D model affects molecular diffusion and dispersion equally it is represented by lowering the effective molecular diffusivity and dispersivity equally in the layered region of the 1-D version instead.Diffusivities are tuned such that the CO 2 pro-Figure 6. Simulated and observed CO 2 and CH 4 concentrations in the firn at WAIS Divide.The model is initialized with the recorded atmospheric trace gas concentrations in 1800 CE at all depths and is forced at the surface with histories of atmospheric CO 2 and CH 4 concentrations (Etheridge et al., 1996(Etheridge et al., , 1998;;Keeling et al., 2005;Buizert et al., 2012;Dlugokencky et al., 2016a, b).Markers indicate observed CO 2 (diamonds) and CH 4 (crosses) concentrations (Battle et al., 2011).Based on high CO 2 and CH 4 , two samples at ∼ 15 and ∼ 50 m depth were likely compromised by modern air during analysis and are thus ignored in the curve fit.Differences in the CO 2 and CH 4 profiles between the 1-D model and the 2-D model with or without barometric pumping are not visible at the resolution of this figure but are illustrated in Fig. S12.
files are (nearly) identical.The small remaining deviations in CO 2 and CH 4 concentrations between model versions (< ±0.4 ppm and < ±9 ppb, respectively) are illustrated in Fig. S12.

δ 15 N and thermal fractionation
In all four model setups, the seasonal cycle of temperature dominates the shape of the δ 15 N profiles in the top ∼ 30 m (Fig. 7).Warm summer temperatures drive the migration of heavy isotopes into the colder firn below and produce a δ 15 N summer peak just below the surface (Severinghaus et al., 2001).In contrast, a minimum of δ 15 N occurred in this region during the previous winter season when the thermal gradients were reversed.The remnants of this winter minimum are still visible in the gas record as anomalously low δ 15 N values below the summer peak.The differences between observations and simulated δ 15 N values in the top of the firn column are likely caused by the idealized representation of our seasonal cycle in the model.

Impact of reduced-permeability layers
In the layered 2-D model without barometric pumping, the simulated δ 15 N values are close to observations at the top of the LIZ but continue to increase with depth (Fig. 7).This is contradicted by observational data even when the unusually low δ 15 N values at the COD and below are not taken into consideration (near the COD, firn air sampling becomes more difficult in the field and the potential for fractionation during sampling is increased).A closer inspection of the LIZ in Figs.7 and 8 reveals a zigzag pattern in the δ 15 N profile where impermeable layers are present.Simulated isotope ratios are higher just above horizontal layers, where heavy isotopes can accumulate, and are anomalously low below layers where heavy isotopes are more readily removed than supplied by gravitational settling.Gravitational settling through gaps in the layers sets up small horizontal concentration gradients that drive horizontal Fickian diffusion.Layering increases the effective travel path length for molecules and reduces the effective vertical diffusivity by increasing the tortuosity.However, layering alone appears to be insufficient to prevent gravitational settling completely, with continued gravitational enrichment being observed in the LIZ in this model version.
Is the impact of layers on the firn trace gas profile larger for isotopes such as δ 15 N than for anthropogenic tracers such as CO 2 , CFCs, or CH 4 ?All three gases have experienced large atmospheric variability in the industrial era.Therefore, the migration of these gases into the firn is dominated by vertical and horizontal Fickian diffusion.For δ 15 N, in contrast, the effect of gravity is critical for transport.To answer the above question, we compare output from the layered 2-D model to the 1-D model without layers.We find that δ 15 N values in the 1-D and 2-D model without barometric pumping are almost identical (Fig. 7).Layering in the 2-D model increases the effective transport distance for CO 2 just as much as for δ 15 N and there is no disproportional impact of layering on gravitationally fractionated isotope ratios.Differences in explicitly represented tortuosity are automatically compensated in the 1-D model during tuning to the same CO 2 profile by reducing molecular diffusivities.Therefore, we conclude that layering alone does not simultaneously explain the observed CO 2 and δ 15 N patterns.

Barometric pumping and the emergence of the LIZ
δ 15 N values simulated by the 1-D and 2-D model with barometric pumping are lower in the LIZ than in both model versions without barometric pumping (Fig. 7).Accounting for barometric pumping improves the agreement with observations throughout the LIZ.However, the reduction of gravitational fractionation is substantially stronger when layers are present.Only when both layering and barometric pumping are accounted for in the model does the δ 15 N profile correctly indicate no further gravitational enrichment in the LIZ and matches observations more closely.Dispersive mixing is nearly independent of molecular mass and does not lead to gravitational fractionation, but rather acts to elimiwww.the-cryosphere.net/12/2021/2018/The Cryosphere, 12, 2021-2037, 2018 nate the concentration gradients associated with gravitational settling.Although barometric pumping velocities are largest near the surface (Fig. 3c), significant dispersive mixing is generally limited to the LIZ because dispersivity of firn is inversely related to the open porosity in the parameterization of Buizert and Severinghaus (2016) and dispersion is overwhelmed by molecular diffusion in the DZ.Furthermore, molecular diffusivities drop rapidly in the LIZ in the model (Fig. 5).Because dispersion provides an additional transport mechanism for trace species, even less molecular diffusion is needed to match observed CO 2 and CH 4 concentrations in the LIZ when barometric pumping is active.Layering amplifies the importance of barometric pumping because gravitational fractionation between annual layers is restricted to the small gaps in the LIZ (Fig. 8).Narrow pathways amplify barometric pumping flow velocities and thus dispersive mixing in these regions (Fig. 4), thus overwhelming the influence of gravitational fractionation more readily than in the 1-D model.This effect is responsible for the larger differences between the δ 15 N profiles obtained from the two model versions with barometric pumping in Fig. 7.The strength of dispersive mixing in our layered 2-D model is physically motivated; thus, barometric pumping and layering together lead to a more natural emergence of the δ 15 N-defined LIZ in the 2-D model.

Surface mixed zone depth
We estimate the depth of the SMZ at WAIS Divide to be ∼ 2.8 m.Multiple different procedures have been used to estimate SMZ thickness in the past, many of which rely on δ 15 N data in the deep firn near the LIZ (Battle et al., 2011). The

Law Dome DSSW20K
At Law Dome DSSW20K, the firn thickness is ∼ 20 m less than at WAIS Divide.Accumulation rates are comparable, but annual-mean temperatures are ∼ 10 K warmer.The SMZ is slightly deeper and barometric pumping is stronger at Law Dome, yielding more near-surface and dispersive mixing.
Constraining the SMZ depth at DSSW20K is more difficult because fewer δ 15 N measurements are available for this site, and their associated uncertainty is, at ±15 per meg, much larger than at the more recently sampled WAIS Divide site (Trudinger et al., 2002).Molecular diffusion generally takes a less important role at DSSW20K and molecular diffusivities obtained by tuning are about half or less than those at WAIS Divide for most of the firn column.Thermal fractionation has a weaker impact on the isotope record near the surface at Law Dome due to the smaller amplitude of the seasonal cycle and stronger near-surface mixing compared to WAIS Divide.Figures of molecular diffusivity, advection velocities, and other firn properties at DSSW20K are provided in the Supplement.
Simultaneously matching the δ 15 N, CO 2 , and CH 4 profile at Law Dome DSSW20K has proven difficult in the past (Trudinger et al., 2002;Buizert and Severinghaus, 2016).Simulated δ 15 N in the LIZ is typically substantially higher than in observations.Buizert and Severinghaus (2016) suggested that barometric pumping in the deep firn may be able to reconcile this contradiction.However, the mixing obtained from theoretical predictions was insufficient to achieve a satisfactory fit.Buizert and Severinghaus (2016) hypothesized that firn layering may play a critical role in amplifying the impact of barometric pumping.The authors used an idealized eddy and molecular diffusivity profile in the deep firn to simulate the effect of layers on firn air transport.Using these diffusivity profiles, they were able to obtain good agreement with observed δ 15 N, CH 4 , and 14 CO 2 .Our 2-D model includes an explicit representation of layering and is subject  (Trudinger et al., 2002).(Trudinger et al., 2002(Trudinger et al., , 2013)).The dashed black line represents the equilibrium solution for pure gravitational settling (δ grav ).The horizontal blue line marks the depth where vertical diffusivity reaches zero.The inset shows a magnification of the lock-in zone.
to similar physical constraints on barometric pumping as the 1-D model of Buizert and Severinghaus (2016)  patterns of both profiles are reproduced correctly (Fig. 9).But the disagreement between modeled and observed δ 15 N in the deep firn remains despite barometric pumping producing substantial non-fractionating dispersive mixing in the region (Fig. 10).Simulated δ 15 N values diverge from observations at ∼ 38 m, where gravitational enrichment seems to stop in observations but continues in the model.In contrast, the LIZ, as defined by CO 2 and CH 4 , only starts at roughly 43 m depth.Such an early onset of dispersive mixing is not supported by the dispersivity parameterization.However, only the longitudinal-to-flow mixing in the vertical direction at Siple Station was used to develop the firn dispersivity parameterization, and the use of this parameterization may be inappropriate at Law Dome DSSW20K (Buizert and Severinghaus, 2016).Moreover, dispersivity typically differs in the horizontal and vertical as well as the longitudinal-to-flow and the traverse-to-flow directions, an effect that is not accounted for in this study because of a lack of observational evidence to constrain anisotropic dispersivity.

Differential kinetic isotope fractionation
Isotope ratios in firn typically do not reach values as high as predicted from gravitational equilibrium due to the influence of advection and non-fractionating dispersive mixing . is defined as the (typically negative) difference between any mass-normalized isotope ratio and δ 15 N (see text).Subscripts of one or two element names identify ratios as isotope or elemental ratios, respectively.The dashed black line highlights where molecular diffusivity in the model reaches zero.(Trudinger et al., 1997;Kawamura et al., 2013;Buizert and Severinghaus, 2016).Advection and mass-independent mixing transport less-fractionated air down in the firn column and act to counterbalance the enrichment of heavy isotopes by gravitational fractionation.As a result, all isotope ratios fall below the gravitational settling line δ grav (Fig. 11) but the magnitude of the deviation depends on the specific isotope pair.
The magnitude of disequilibrium of different isotope and elemental ratios is quantified here by defining the (massnormalized) kinetic fractionation relative to δ 15 N (following Kawamura et al., 2013) as where m x/y is the mass difference (in mol kg −1 ) of isotopes x and y.This definition is similar to the 86 Kr excess terminology introduced by Buizert and Severinghaus (2016) but is given in the more precise ln(q) notation and uses δ 15 N as reference instead of δ 40 Ar / 36 Ar.To calculate x/y , isotope ratios must have been previously corrected for the influence of thermal fractionation by removing temperature effects with a suitable firn air transport model (Fig. 12).The degree of disequilibrium, represented by , is controlled by differential kinetic isotope fractionation.Heavy, slow-diffusing isotopes approach gravitational equilibrium more slowly than lighter, faster-diffusing isotopes.Therefore, slow-diffusing isotopes experience larger kinetic frac- The Cryosphere, 12, 2021Cryosphere, 12, -2037Cryosphere, 12, , 2018 www.the-cryosphere.net/12/2021/2018/tionation (i.e., deviations from gravitational equilibrium) in regimes with non-zero advection or dispersion.Consequently, is more negative for heavier, slower-diffusing isotopes.On its own, this rule of thumb cannot fully explain the pattern of ratios containing two different elements, such as 132 Xe / 28 N 2 .The magnitude of disequilibrium in such mixed-element ratios is further discussed in Appendix B.
In the DZ, decreases almost linearly with depth, while within the SMZ and LIZ changes much more rapidly.Where molecular diffusivity is zero, remains constant.This pattern is explained by the relative importance of advection and dispersive mixing compared to molecular diffusion in different regions of the firn column.The Péclet number (Pe) traditionally quantifies the ratio of the advective to the diffusive transport and is here defined as the ratio of the diffusive to the advective timescale (τ D m and τ adv ).We add the timescale of dispersive mixing (τ D e ) to the numerator because the effect of advection and dispersive mixing on the isotope profiles is very similar although the physics differ (Kawamura et al., 2013): where L = 1 m is the characteristic length scale of firn air transport and D m , W , and D e are characteristic values of the molecular diffusivity, the time mean vertical advection velocity, and the vertical dispersive or convective mixing at that depth.
This modified Péclet number varies in the model by several orders of magnitude through the firn column at WAIS Divide with peak values in the SMZ and the deep firn (Fig. 13).High modified Péclet numbers in the SMZ are caused primarily by large D e values, and high modified Péclet numbers in the LIZ are mostly the result of low molecular diffusivities.Kawamura et al. (2013) showed analytically that relative kinetic isotope fractionation depends on the ratio of eddy diffusivity to molecular diffusivity, but the role of advection was neglected due to near-zero accumulation rates at the Megadunes site.The absolute difference in kinetic isotope fractionation (i.e., ) should be greatest when the product of the modified Péclet numbers of both isotopes is near one.In line with these theoretical predictions, we observe almost no further isotopic enrichment of δ 15 N in the LIZ when barometric pumping is included in the model and Pe 1 (Figs.7, 13).The largest changes of occur in the 2-D model when the product of the modified Péclet numbers is within approximately 1-2 orders of magnitude of unity.This region is illustrated by the vertical gray bar in Fig. 13, which contains the SMZ as well as the beginning of the LIZ where non-fractionating mixing is of similar magnitude as molecular diffusion.
With active barometric pumping and centimeter-scale layering, the product of the modified Péclet numbers at the bottom of the LIZ becomes so large that stops to decrease entirely in our model.If barometric pumping is neglected instead, the modified Péclet numbers in the layered 2-D model are considerably lower in the LIZ and some gravitational and kinetic fractionation persists (i.e., δ 15 N and continue to change gradually).Therefore, barometric pumping leads to slightly weaker rather than stronger differential kinetic fractionation at the COD of WAIS Divide in the model, in contrast to expectations (Buizert and Severinghaus, 2016).Furthermore, layering and barometric pumping in the model seem to be insufficient to obtain the full ∼ 5-23 per meg per amu range of Kr values measured in the WAIS Divide ice core record (WAIS Divide Project Members, 2015;Bereiter et al., 2018). Instead, other, unresolved (i.e., subgrid-scale) processes may be the reason for the large observed disequilibrium.Establishing a straightforward relationship between disequilibrium and surface pressure variability using firn air models alone may not be possible without more observational data.

Diffusive fractionation
Strong kinetic isotope fractionation can also be observed for trace gases that experience large changes in the atmospheric mixing ratio while the atmospheric isotope ratios remain constant (Trudinger et al., 1997;Buizert et al., 2013).As the concentration of a trace gas increases, isotopologues of the gas migrate into the firn column at different speeds because of small differences in their masses and diffusivities.This results in a relative depletion of the slower-diffusing isotopologue with depth called diffusive fractionation (Trudinger et al., 1997).During periods of abrupt CH 4 change, diffusive fractionation commonly amounts to a notable correction in ice core studies (Trudinger et al., 1997;Buizert et al., 2013).Diffusive fractionation of δ 13 C-CH 4 is strong and poorly constrained by models, to the degree that it prohibits the reliable atmospheric reconstruction of this parameter from firn air measurements (Sapart et al., 2013).Since diffusive fractionation is a type of kinetic fractionation, it can be tested in our model.We assume a constant atmospheric 13 C / 12 C isotope ratio of 1.1147302 % for CO 2 (i.e., δ 13 C-CO 2 = −8 ‰) and 1.0709052 % for CH 4 (i.e., δ 13 C-CH 4 = −47 ‰).Thermal fractionation and gravitational settling are neglected to isolate the impact of the atmospheric mixing ratio change.The model including barometric pumping calculates δ 13 C-CO 2 and δ 13 C-CH 4 values depleted by up to ∼ 0.2 and ∼ 2 ‰ relative to the atmosphere in the WAIS Divide LIZ at the time of firn air sampling (Fig. 14).Without barometric pumping, delta values are notably higher because molecular diffusion is stronger, and the dispersive mixing no longer smooths out the profile in the deep firn.

Predicting disequilibrium
Past mean ocean temperature can be estimated from the noble gas concentrations in ice core bubbles (Headly and Severinghaus, 2007;Ritz et al., 2011;Bereiter et al., 2018).On glacial-interglacial timescales, atmospheric concentrations of noble gases are primarily controlled by gas dissolution in the ocean.Because the temperature sensitivity of solubility is different for each gas, measurements of noble gas ratios in ice cores can be used to obtain a signal of integrated ocean temperature.However, as for any gas, the trace gas concentrations in bubbles must first be corrected for alterations of the atmospheric signal in the firn.In a recently published deglacial mean ocean temperature reconstruction, the WAIS Divide noble gas ice core record was corrected for gravitational fractionation and thermal fractionation using δ 40 Ar / 36 Ar measurements and a firn temperature gradient estimate (Bereiter et al., 2018).The authors further noted that different degrees of deviation from gravitational equilibrium (i.e., disequilibrium) can bias the gravitational fractionation correction applied to the raw noble gas record, which may lead to a cold bias of ∼ 0.3 • C in Holocene and Last Glacial Maximum temperatures.Disequilibrium effects of −287.5 per meg for Kr / N 2 , −833.3 per meg for Xe / N 2 , and −545.7 per meg for Xe / Kr simulated by our model correspond to absolute temperature biases of 0.33, 0.41, and 0.45 • C, respectively, following the method of Bereiter et al. (2018).Because the magnitude of disequilibrium depends on firn properties and accumulation rate, glacial-interglacial changes in environmental boundary conditions may also affect the magnitude of disequilibrium in firn and thus the size of the relative deglacial temperature change estimated from noble gases.
In an attempt to compensate for disequilibrium effects and gravitational settling at the same time, it has been suggested that the elemental ratios Kr / N 2 and Xe / N 2 in bubbles should be corrected by subtracting krypton or xenon isotope ratios, respectively (Headly, 2008).This would assume that krypton and xenon isotopes may be influenced similarly by the processes responsible for creating disequilibrium in Kr / N 2 and Xe / N 2 .Therefore, this approach may compensate for disequilibrium effects and gravitational settling simultaneously, but it has been untested in models so far.values simulated here allow us to test this method quantitatively.We use a linear fit to predict Kr / N 2 from Kr .The linear fit yields good agreement with the modeled Kr over the firn column (R 2 > 0.998), indicating that the scaling between values is nearly independent of depth.We find that (mass-normalized) Kr / N 2 should be approximately 75 % of (mass-normalized) Kr at the WAIS Divide site.Scaling relationships for other isotope and element pairs are shown in Table 1 and are equally robust.Moreover, our results show that the source of disequilibrium is irrelevant for the correction for the macroscopic processes represented in our model.Advection and convective or dispersive mixing all show the same scaling relationships for .At Law Dome DSSW20K, the calculated ratio of Kr / N 2 and Kr is at 75.9 % almost identical to the result at WAIS Divide.Sensitivity tests with the 1-D analytical model presented in Appendix A demonstrate that the disequilibrium scaling relationship between Kr isotopes and Kr / N 2 is robust to within ±5 % over a wide parameter range of molecular diffusivity, eddy diffusivity, and advection velocity.Uncertainties become largest in the extreme case when Ar , the lowest simulated value, is used to predict Xe / Kr , the highest simulated value, but they never exceed ±25 %.
This suggests that the same scaling relationship between Kr / N 2 and Kr may be assumed to hold for any ice core site without introducing large biases.Kr and Xe from combined measurements of δ 15 N, δ 86 Kr / 82 Kr, and δ 136 Xe / 129 Xe in ice cores could be used to predict the disequilibrium effects on noble-gas elemental ratios (i.e., Kr / N 2 , Xe / N 2 , and Xe / Kr ) and allows us to make a gasspecific gravitational correction.Although predicted Kr values at WAIS Divide are close to the current analytical uncertainty of the 86 Kr / 82 Kr measurement, correcting for kinetic fractionation and disequilibrium will become advisable with future improvements in precision and may improve mean ocean temperature reconstructions.

Conclusions
We developed a two-dimensional firn air transport model that explicitly represents tortuosity in the firn column through migrating layers of reduced permeability.The idealized representation of firn layering is physically motivated and may illustrate the impact of firn density anomalies (i.e., summer vs. winter firn or wind crusts) on gas transport.The model also accounts for thermal fractionation, a surface mixed zone, and surface pressure-forced barometric pumping.Dispersive mixing as a result of barometric pumping is constrained in the model by previously published parameterizations and not subject to tuning.Simulations of the δ 15 N profile at WAIS Divide show that extensive horizontal diffusion through the tortuous firn structure is required by the discontinuous layers.This limits the effective vertical diffusion of gases at depth.However, layering alone does not sufficiently reduce gravitational enrichment of isotopes in the deep firn.Similarly, the effect of barometric pumping alone is insufficient to obtain agreement with observations.The combination of barometric pumping with layering, in contrast, leads to amplified dispersive mixing.This is due to high velocity focusing in layer openings and leads to a more natural emergence of a lock-in zone in the model.
Previous studies have shown that downward advection, convective mixing, and dispersive mixing all hinder trace gases in reaching the isotope ratios expected from gravitational settling (e.g., Severinghaus et al., 2010;Kawamura et al., 2013;Buizert and Severinghaus, 2016).Kinetic fractionation is strongest for slow-diffusing gases and increases with firn column depth.Our numerical experiments show that barometric pumping leads to increased isotopic disequilibrium in the firn column.However, our simulations fail to account for the full range of 86 Kr excess observed in the WAIS Divide core, as well as for the relatively weak δ 15 N enrichment seen at DSSW20K, suggesting that these effects are not caused by the presence of layering (as previously suggested) and that their origin must be sought elsewhere.We further find robust scaling relationships between the magnitude of disequilibrium in different noble gas isotope and elemental ratios.Our results suggest that, to first order, these scaling relationships are independent of depth in the firn column and independent of the cause of disequilibrium.Thus, a correction that accounts for differential kinetic fractionation may be applied to observed noble gas ratios in the reconstruction of mean ocean temperature.Analysis of this relationship reveals that disequilibrium should most strongly affect ratios of two heavy isotopes, such as 132 Xe / 84 Kr, because heavy elements diffuse more slowly than N 2 (i.e., x/28 0) and the mass weighting factor is larger in the first than in the second term (i.e., m x −m 28 m x −m y m y −m 28 m x −m y ).Although this equation can theoretically predict of any isotope ratio from of the two isotopes x and y relative to 28 N 2 (i.e., x/28 and y/28 ), in practice, this approach will not allow correcting for differential kinetic isotope fractionation.x/28 cannot be measured directly and the atmospheric ratio of the noble gas x to nitrogen is not constant over long timescales.Thus, x/28 will not only be affected by disequilibrium but will also be influenced by atmospheric variability resulting from gas specific solubility differences (i.e., precisely the mean ocean temperature signals we attempt to reconstruct).Instead we suggest that the scaling relationships provided in Sect.4.3 can be used to predict the of noble gas elemental ratios.

2022B.
Figure 1.Layering of firn photographed in a surface pit at WAIS Divide.Image courtesy of Anaïs Orsi.

Figure 2 .
Figure 2. Schematic depiction of a typical isotope profile.The surface mixed zone (SMZ), the diffusive zone (DZ), the lock-in zone (LIZ), and the ice below are indicated by shading.Further indicated are the lock-in depth (LID) and the close-off depth (COD) (see text).

Figure 4 .
Figure 4.The different components of the velocity field: (a) firn velocity, (b) the velocity of air return flow to the atmosphere due to pore compression, and (c) net firn air advection velocity (a linear combination of the fields in panels a and b).Because of its alternating direction, barometric pumping yields no net flow but instantaneous flow field patterns look similar to panel (b).Black arrows indicate the downward advection of layers at the firn velocity.
Figure 5.The CO 2 diffusivity profile at WAIS Divide.(a) Horizontally averaged, vertical and horizontal diffusivity in the model with and without barometric pumping.(b) Map of diffusivity in the 2-D model without barometric pumping.Only every third layer present in the model is shown here for clarity.

Figure 7 .
Figure 7. Horizontally averaged δ 15 N at WAIS Divide.Model output is shown from four different versions of the 2-D model (see text).Black circles with error bars indicate the observed firn δ 15 N (Battle et al., 2011).The dashed black line represents the equilibrium solution for pure gravitational settling (δ grav ).The horizontal blue line marks the depth where vertical diffusivity reaches zero.The inset shows a magnification of the lock-in zone.

Figure 8 .
Figure 8. Simulated δ 15 N in a section of the lock-in zone at WAIS Divide from the 2-D model not including barometric pumping.Impermeable horizontal layers are shown in red.The size of the openings in the layers shrink with increasing depth.
Cryosphere, 12, 2021-2037, 2018 www.the-cryosphere.net/12/2021/2018/B. Birner et al.: Layering and barometric pumping in firn 2029However, if the deep firn is affected by dispersive mixing due to barometric pumping, these estimates may be falsely attributing some fraction of the dispersive mixing in the deep firn to the SMZ.To address this problem, we follow the method ofSeveringhaus et al. (2001) in calculating SMZ thickness.This approach compares the depth where δ 15 N reaches a certain value in two different model configurations with zero and non-zero values of D SMZ .Thermal effects are neglected.The first setup is the 2-D model with barometric pumping as presented above but the dispersivity is set to zero everywhere without retuning the model.The SMZ thickness is calculated from the depth difference between this model run and a second model run where barometric pumping is deactivated and D SMZ = 0. Thus only advection and gravitational fractionation shape the profile of δ 15 N. Our depth estimate of 2.8 m is within the 1.4-5.2m range published previously(Battle et al., 2011).

Figure 9 .
Figure 9. Simulated CO 2 and CH 4 concentrations in the firn at Law Dome DSSW20K.The model is forced with histories of atmospheric CO 2 and CH 4 concentrations from 1800 to 1998 CE (the date of sampling).Markers indicate observed CO 2 (diamonds) and CH 4 (crosses) concentrations(Trudinger et al., 2002).

Figure 10 .
Figure 10.Horizontally averaged δ 15 N at Law Dome DSSW20K.The solid lines show the results of the 2-D model (with layers) for the cases with (blue) and without (red) barometric pumping.Black circles with error bars indicate the observed firn δ 15 N(Trudinger et al., 2002(Trudinger et al., , 2013)).The dashed black line represents the equilibrium solution for pure gravitational settling (δ grav ).The horizontal blue line marks the depth where vertical diffusivity reaches zero.The inset shows a magnification of the lock-in zone.

Figure 11 .
Figure 11.Horizontally averaged isotope ratios at WAIS Divide in the 2-D model including barometric pumping and horizontal layers.Isotope ratios are normalized to one atomic mass unit (amu) mass difference (Supplement).The dashed black line represents the equilibrium solution for pure gravitational settling (δ grav ).The horizontal blue line marks the depth where vertical diffusivity reaches zero.Observed δ 15 N are shown as circles with horizontal error bars(Battle et al., 2011).The inset shows a magnification of the lock-in zone (gray patch).

Figure 12 .
Figure12.Differential kinetic isotope fractionation ( ) profiles for different isotope pairs at WAIS Divide.Colored solid and dashed lines show results from the 2-D model with and without barometric pumping, respectively.is defined as the (typically negative) difference between any mass-normalized isotope ratio and δ 15 N (see text).Subscripts of one or two element names identify ratios as isotope or elemental ratios, respectively.The dashed black line highlights where molecular diffusivity in the model reaches zero.

Figure 13 .
Figure 13.The balance of fractionating and non-fractionating mixing at WAIS Divide.Panel (a) illustrates the horizontally averaged modified Péclet number of δ 15 N (Pe 15 , see text).Blue and red lines show results from the 2-D models with and without barometric pumping.The strength of dispersive mixing in the calculations is given by the mean barometric pumping flow velocities at the site.Panel (b) displays the product of the modified Péclet numbers for δ 15 N and 84 Kr / 28 N 2 (Pe 84 ).The region where changes with depth should be greatest (i.e., where the product of the modified Péclet numbers is near 1) is highlighted by a gray bar.Panel (c) provides a magnified 2-D map of this Péclet number product in the LIZ.

Figure 14 .
Figure 14.Diffusive fractionation effect at the time of sampling at WAIS Divide on δ 13 C of (a) CO 2 and (b) CH 4 .Atmospheric mixing ratios of 12 CO 2 , 13 CO 2 , 12 CH 4 , and 13 CH 4 were obtained from atmospheric trace gas histories used to drive the firn air model (see Supplement) and assuming constant atmospheric isotope ratios of −8 and −47 ‰ for δ 13 C-CO 2 and δ 13 C-CH 4 , respectively.Firn air values are presented as the difference from the constant atmospheric isotope ratios.

Table 1 .
The scaling factor in the 2-D model with barometric pumping between different element and isotope ratios from linear regression of value pairs at all depths.R 2 > 0.996 for all relationships.