A new multi-gas constrained model of trace gas non-homogeneous transport in firn : evaluation and behavior at eleven polar sites

Insoluble trace gases are trapped in polar ice at the firn-ice transition, at approximately 50 to 100 m below the surface, depending primarily on the site temperature and snow accumulation. Due to the different time scales for snow accumulation versu s diffusion of gases through the snowpack, age differences between gases and the ice in which they are “trapped” can be large; e.g. several thousand years in central Antarctica (a low snow accumulati on area). Models of trace gas transport 5 in polar firn are used to relate firn air and ice core records of t race gases to their atmospheric history. We propose a new model based on the following contributions. Fir t, the airflow transport model is revised in a poromechanics framework with emphasis on the no n-homogeneous properties and the treatment of gravitational settling. We then derive a non-l inear least square multi-gas optimization scheme to calculate the effective firn diffusivity (automat ic diffusivity tuning). The improvements 10 gained by the multi-gas approach are investigated (up to ten gases for a single site are included in the optimization process). We apply the model to four Arctic (De von Island, NEEM, North GRIP, Summit) and seven Antarctic (DE08, Berkner Island, Siple Dome, Dronning Maud Land, South Pole, Dome C, Vostok) sites and calculate their respective depthdependent diffusivity profiles. Among these different sites, a relationship between an increasin g th ckness of the lock-in zone defined from 15 the isotopic composition of molecular nitrogen in firn air (d enotedδN) and the snow accumulation rate is calculated, in accordance with observations. It is a ssociated with reduced diffusivity depthgradients in deep firn, which decreases gas density depth-gr adients, at high accumulation rate sites. This has implications for the understanding of δN of N2 records in ice cores, in relation with past variations of the snow accumulation rate. Although the snow accumulation rate is clearly a primary 20 control on the thickness of the lock-in zone, our new approac h that allows for the calculation of an estimated lock-in depth may lead to a better constraint on th e age difference between the ice and entrapped gases.


Introduction
Modelling of gas transport in firn has been an active field of research in the past two decades, with physical models developed and refined in parallel with measurement campaigns. Without pretending to be exhaustive, several key results can be mentioned here. The enclosure of air in bubbles associated with firn sinking was quantified by measurements of closed porosity and modelled by Stauffer et al. (1985). Gas transport has first been described as a molecular diffusive process with gravitational correction by Schwander (1989). A mathematical model of the impact of the firn sinking velocity and bubble trapping was introduced simultaneously by Rommelaere et al. (1997) and Trudinger et al. (1997). Rommelaere et al. (1997) also introduced calculation of trace gas concentration in closed bubbles. The impact of thermal fractionation and turbulent transport in the convective layer (eddy flows) on the diffusion phenomenon was described by Severinghaus et al. (2001). These models were primarily focused on a diffusive process driven by a concentration gradient (Fick's law). The impact of permeability on the transport model was investigated by Schwander (1989) and Freitag et al. (2002).
Modelling firn air transport for the purpose of reconstructing a trace gas atmospheric history can be decomposed into three steps. First, a physical transport model (often referred to as the forward model) describes the gas behaviour in the ice lattice, including the impact of medium heterogeneities (depth-dependent porosity, gas trapping in bubbles, localised eddy flows, etc.). Six forward models of trace gas transport in firn have recently been inter-compared using a common set of input parameters . The effective diffusivity appears as a key depth dependent parameter to characterise a given firn. This parameter includes all diffusive phenomena experienced by the gas and is simply referred to as diffusivity in the following discussion.
The second step is to calculate the depth profile of diffusivity from measured gases in the firn and their known atmospheric history. The associated model is defined by the combination of the forward model, gas history, borehole measurements and a diffusivity optimisation method, and is typically termed as the inverse diffusivity model, as it calculates the solution of an inverse problem (see e.g. Rommelaere et al., 1997;Trudinger et al., 2002;Buizert et al., 2012). This step provides an optimised diffusivity profile that is consistent with the firn-air observations and known past atmospheric changes. The fact that it can be significantly different from the diffusivity profile measured on small 4 cm samples supports the findings of Fabre et al. (2000), who concluded that firn diffusivity is more influenced by macro-scale structure than by micro-scale features. A major goal of the paper is to derive a diffusivity profile based on multiple gases with known past histories. Solving the associated nonlinear optimisation problem involves numerous runs of the forward model for each species with trial diffusivity profiles, in order to infer the transport model sensitivity to small diffusivity variations. This iterative process puts constraints of computational efficiency and numerical robustness on the forward model algorithm.
The third and last step is to reconstruct the atmospheric history of gases for which only measurements of firn air are available, possibly correlating the results obtained at different polar sites. Indeed, reconstructing the mixing ratio trends of trace gases prior to their atmospheric measurement period is a major motivation for firn air analysis (see e.g. Montzka et al., 2011, and references therein). Several methods were developed for that purpose, based on Green functions (Rommelaere et al., 1997), Monte-Carlo simulations (Bräunlich et al., 2001), or effective age (Trudinger et al., 2002). These methods were recently applied, for example, to HFC-227ea (1,1,1,2,3,3,3-heptafluoropropane, Laube et al., 2010), mercury (Faïn et al., 2009) and HFC-23 (trifluoromethane, Montzka et al., 2010), respectively. Alternatively, the mixing ratio trends of trace gases can be assessed from their atmospheric budgets and compared to firn data using a forward firn model (e.g. Butler et al., 1999;Martinerie et al., 2009;Montzka et al., 2010).
Our goals in this paper are to analyse the transport phenomena associated with the forward model, in order to build an inverse diffusivity model based on multiple gases. This inverse model is evaluated at eleven polar sites. The modelling constraints and assumptions are motivated by the perspective of reconstructing the atmospheric histories of gases using multi-sites inversion methods.
The forward model includes the diffusion process and gravitational settling in the treatment of firn gas transport in a poromechanics framework (interconnected network composed of the ice lattice, the gas connected to the surface and the gas trapped in bubbles) and is described in Sect. 2. The problem formulation and optimisation issues associated with the inverse diffusivity model are then considered in Sect. 3, along with an evaluation of the efficiency of the method and the impact of key parameters. The diffusive behaviour of eleven modelled polar sites are finally compared in Sect. 4, showing a dependency of the deep firn physics on the snow accumulation rate.

Forward model and physical concepts
Mathematical modelling of fluid transport implies to simultaneously solve the time-varying partial differential equations (e.g., Euler or Navier-Stokes equations, see Anderson, 1991, for more details) representing: (1) the conservation of mass (continuity), which relates the time-evolution of the fluid mass or concentration to the flux induced by the rate of motion of the fluid particles (referred to as the fluid velocity and expressed in metres per year in our model); (2) the conservation of momentum (derived from Newton's second law), which describes the effect of internal (e.g., induced by a local pressure gradient) or external (e.g., gravity) forces on the rate of change of the fluid velocity; (3) the conservation of energy, which provides, for example, the temperature evolution. Models for trace gas transport in firn (involving gas and ice) focus on the continuity equation to compute the mixing ratio of trace gases in air. Algebraic formulas for the conservation of momentum (e.g., from Fick's law or using hydrostatic distributions) provide approximations of the average fluid velocity. We do not use the conservation of energy, considering that the firn is isothermal (at mean annual temperature). Based on this framework, our forward model and physical assumptions are, thus, focused on concentration and velocity computation for ice and gases.
The forward model developments are also motivated by requirements of the optimisation procedure for calculating firn diffusivity (Sect. 3). First, as the optimisation procedure of the inverse diffusivity model iterates on (possibly nonphysical) diffusivity profiles, the forward model has to be robust enough (e.g., using smooth functions and avoiding singularities) to provide relevant results with such trial profiles. Second, the nonlinear effects have to be avoided (if possible) to prevent the existence of multiple solutions to the optimisation problem.

From the firn network to trace gas transport
The conservation of mass, both for ice and gases, can be described by considering the firn as composed of a sinking solid ice structure (i.e., increasing in depth due to continued accumulation of snow at the surface) and of a transported fluid (air and trace gases). The space available for the fluid is composed of the open pores (having a free path to the surface) and the closed pores (sinking with the firn), which constitute a network of two spaces interconnected through the rate of fluid mass transferred from one space to the other. The dynamics of such a system is inferred from mass conservation in a porous media, expressed in a 1-D (vertical) Eulerian frame (attached to the snow surface). The 1-D representation implies a layer-averaged description of the transport phenomena (i.e., the firn has horizontally-homogeneous properties). The major depth-variation of the transport parameters (vertical variation due to e.g., firn densification processes or seasonality) can be represented by such description. It results from variations of the grain size and pore space, which were observed to occur mostly in the vertical direction in polar firn by e.g., Hörhold et al. (2009) using X-ray microtomography. The main limitation of the 1-D assumption is the lack of representation of the effects of cracks or firn layering characteristics (horizontal layers width, extent, physical properties).
The temporal changes in physical properties of the ice lattice, fluids (gases) in open and closed pores are respectively calculated according to Coussy (2003) from the continuity equations: where ρ ice is the ice density (kg per m 3 of void space volume), ρ o gas the gas concentration (mol per m 3 of void space volume) in open pores, ρ c gas the gas concentration in closed pores, z the depth, the total porosity (pores volume/firn layer volume), f the open porosity (open pores volume/firn layer volume), v the firn sinking velocity (rate of increase in depth, in m/year), w gas the relative gas velocity with respect to the firn and r o→c the rate of gas mass transport from the open to the closed pores. The "gas" subscript generically denotes the gas considered (air or trace gas in our case). ∂/∂t and ∂/∂z denote the partial derivatives with respect to time and depth, respectively. The main notations are summarised in Table 1 of the Supplement. A radioactive decay with rate λ can also be considered in the case where a radiogenic substance (e.g., 14 CO 2 in this study) is the tracer. Equations (1a) and (1c) describing the ice lattice and gases mixing ratio in closed pores are solved as in Rommelaere et al. (1997). Our focus is on the gas concentration in open pores (Eq. 1b) and the role of diffusivity in trace gas transport.
In order to calculate the firn sinking velocity and open porosity from the available data, the following assumptions are necessary: (1) the temperature and pressure variations have a negligible effect on the density of pure ice (ρ ice is constant), (2) the variations of the snow accumulation rate have a negligible impact on the ice lattice sinking and the firn densification is assumed to be constant (Eq. (1a) is considered at steady-state), (3) the firn bulk density profile ρ firn (z) is smooth and monotonically increasing with depth (the layering is neglected), (4) the firn thickness and surface elevation remain constant. The firn sinking velocity v(z) = a accu /[ρ ice (1 − (z))] is inferred from the surface accumulation rate of snow a accu (kg m −2 yr −1 ). The total porosity is calculated from the firn density ρ firn (z) and ρ ice , with the temperature-dependent ice density proposed by Schwander et al. (1997), as: where T is the site temperature. In the computation proposed by Goujon et al. (2003), the open porosity f is set by the mean close-off density ρ co (density at which f/ = 0.37).
Here this calculation is slightly modified as in Buizert et al. (2012): ρ co is obtained by specifying the full close-off depth z F (depth at which all the pores are closed, i.e., f (z F ) = 0) independently for each site and solving: E. Witrant et al.: Multi-gas transport in firn at eleven polar sites ρ co = ρ ice − (ρ ice − ρ firn (z F )) × (1/0.37) 1/7.6 (4) The rate of gas mass exchange between the open and closed networks is modelled as a loss rate term r o→c gas = τ (z)ρ o gas , where the trapping rate τ (z) = −v ∂[f/ ]/∂z takes into account the firn sinking and open versus closed porosity, as proposed by Rommelaere et al. (1997). The selective permeation of gases (depending on their effective molecular diameters and occurring during the bubble closeoff) proposed by Severinghaus and Battle (2006) is not accounted for. However, this process is only significant for atoms and molecules of small diameter, which are not considered in our study.

An advective-diffusive model of trace gas transport in open pores
Considering the transport of trace gases in open pores, our firn modelling approach allows separating the impact of gas trapping and advection with firn sinking from the impact of gas particles motion in the volume of open pores. More precisely, denoting a specific trace gas with the subscript α, its concentration variations in open pores is determined from Eq. (1b) as: with ρ o α (0, t) = ρ atm α (t), the trace gas concentration in the overlying atmosphere. The influence of the rate of motion of gas particles on the evolution of ρ o α is included with the trace gas velocity with respect to the firn w α (z, t). A direct approach to solving for Eq. (7) would require to calculate the time-variation of w α (see e.g. Coussy, 2003). Instead, a simplified quasi steady-state description of w α is proposed below, in order to limit the number of differential equations.
Considering our target components as trace species in air and that the gases do not affect each other (e.g., no chemical reaction), the classical advective-diffusive model (see e.g. Scanlon et al., 2002) distinguishes the advection of the trace gas with air from its molecular diffusion within the air column. The average velocity of a trace gas is, thus, decomposed as: where w air is the advection-driven velocity (with air) and w α = (w α − w air ) is the rate of trace gas motion in the air column induced by molecular diffusion (described by Fick's law). This model also supposes no interaction between the diffusive processes for each gas.
The air velocity can be calculated at steady-state (no timevariation, ∂ρ o air /∂t = 0 in Eq. 1b) by solving the boundary value problem: The steady-state air concentrationρ o air is approximated by the hydrostatic equilibriumρ o air (z, t) = ρ atm air (t)e M air gz/RT , where ρ atm air is the air concentration in the overlying atmosphere, R the ideal gas constant, g the acceleration due to gravity and M air the molar mass of air. The air velocity w air is induced by the presence of firn sinking (first term in Eq. 9) and bubble trapping (last term, involving also a possible radioactive decay). A similar approach was proposed by Rommelaere et al. (1997).
The relative motion of a trace gas α with respect to air is affected by molecular diffusion according to Fick's first law: where χ air and χ α are the mole fractions of air and gas α, respectively, and D α (z) is the effective molecular diffusivity of gas α. Air being the ambient gas (χ air ≈ 1 and χ α ≈ M air ρ α /M α ρ air ) and considered at steady-state, Fick's law can be written in terms of concentrations as: The relative motion of a gas in the air column associated with molecular diffusion (left term in the equation) is, thus, driven by the gas and air concentration gradients, weighted by the effective diffusivity (right term). An equilibrium set by Fick's law (no diffusive flux) computed as in Eq. (11) implies a constant mixing ratio ρ α /ρ air in the firn column (i.e., w α = 0 in that case). The gravitational settling of a gas in the air column is included by considering instead the flux induced by the difference of the air and gas velocities from their velocities at the steady-state equilibrium. The resulting trace gas velocity is, thus: where the steady-state gas velocity in air w α is directly obtained by using the gas concentration at steady-stateρ o α instead of ρ o α in Eq. (11). Approximating the gas steadystate behaviour by its hydrostatic equilibrium asρ o α (z, t) = ρ atm α (t)e M α gz/RT , where M α is the trace gas molar mass, the advective-diffusive behaviour of the trace gas is modelled as: Considering the continuity Eq. (7) with the trace gas velocity decomposed as in Eq. (13) and the air velocity obtained from Eq. (9), the trace gas concentration in open pores becomes: where the first term in the first equation sets the time variation, the second term is the advection induced by the firn sinking and air velocity, the third term is a mass loss rate due to gas trapping in bubbles (and radioactive decay, where relevant), and the right term is the diffusive transport and gravitational settling. The second and third equations (boundary conditions), respectively, set the solution of the transport equation according to the gas atmospheric concentration and a no-flux condition at the firn bottom z F . The diffusive component in Eq. (14) is equivalent to the diffusion flux induced by gravity proposed by Schwander et al. (1993) which is included in firn transport models (e.g. Rommelaere et al., 1997;Trudinger et al., 1997), as discussed in the Supplement. Schwander et al. (1993) described this term as the influence of gravity on the diffusion model while Rommelaere et al. (1997) referred to this effect as a gravitational correction, but the specific inclusion of gravity using Fick's law was not discussed. The above analysis shows how gravitational settling can be obtained from mass conservation and Fick's law by defining an hydrostatic equilibrium and considering the diffusive flux induced by the deviation from this equilibrium (quasi-steady state description). Using gas velocities (for air and trace gas) as a major variable to describe the gas transport in air allowed us to include both advection and diffusion in this quasi-steady state framework. The respective importance of the different terms in Eq. (13) is evaluated with experimental results in Sect. 2.4.

Turbulent flows and eddy diffusivity
The molecular diffusivity is related to the gas diffusion coefficient in free air D α,air as D α = f D α,air /ν, where ν(z) is the tortuosity (e.g. Coussy, 2003). The molecular diffusivities of two different gases α and β are related by D α /D β = D α,air /D β,air . Thus, molecular diffusivities of trace gases can be simply scaled to the values for a reference species (generally CO 2 ). The flow in the upper layers of the firn may experience fast transients due to the suction effect of the wind, pressure variations, seasonal temperature gradients and the higher porosity. This behaviour is mainly turbulent and the gas dynamics in the upper region can be modelled by introducing an eddy component D eddy (z) in the diffusivity description, determined by the nature of the flow field (e.g., see Byron Bird et al., 2007, for the analogy with Fick's first law and inclusion in the continuity equation). An eddy term characterised by an exponential decay with depth which amplifies (doubles or triples) the effective diffusivity in comparison with the molecular diffusivity to model a convective zone was proposed by Severinghaus et al. (2001). The corresponding eddy flux was introduced in a trace gas transport model in Severinghaus and Battle (2006).
As an alternative approach, we set the maximum admissible molecular CO 2 diffusivity to its value in free air D CO 2 ,air . The effective eddy diffusivity depth z eddy is then defined as the depth above which D α ≥ r α c f D CO 2 ,air , where r α denotes the relative diffusivity of gas α with respect to that of CO 2 in air. The correction factor c f represents the influence of open porosity on the maximum molecular diffusivity in firn and is chosen as an approximation of the surface value of f . For the purpose of our model, the important issue is that D eddy does not depend on the gas considered and should be much larger than the molecular diffusivity (Anderson, 1991). This phenomenon can be modelled by introducing an extra term in the diffusivity that does not depend on the specific gas considered and the overall diffusivity is defined as: where D CO 2 is the CO 2 diffusivity profile in the firn. In this frame, the quantity minimised by the inverse algorithm for diffusivity calculation (described in Sect. 3) is D eddy above z eddy and D CO 2 below z eddy . Several formulations of the eddy region are tested in the next section.

Multi-gas transport at NEEM and Vostok
Our forward model is evaluated at two sites: NEEM, which is the best constrained in terms of number of available trace gases for diffusivity calculation , and Vostok, at which a deep convective zone is observed (Bender et al., 1994). The datasets associated with both sites are described in the Supplement. The choice of these two sites is also motivated by their differences in terms of surface temperature and accumulation rate, thus, providing a large spectrum of surface climatic conditions for testing the model behaviour. The model evaluation involves the calculation of an optimal diffusivity profile (described in Sect. 3). Two types of quality indicators are used here: (1) the root-mean-square deviation (RMSD defined as in Eq. 19) between model results and firn data for the species with varying atmospheric trends used in the diffusivity calculation (reference gases, excluding isotopic indicators as discussed in Sect. 3.3), and (2) the comparison of model results with isotopic indicators of firn air transport δ 15 N, δ 40 Ar and δ 86 Kr 1 . Indeed, as these gas species have a constant atmospheric isotopic composition the repartition of isotopes in the firn is imposed by fractionation related to the mass difference between the major Table 1. Root-mean-square deviation between model results and firn data (RMSD, defined in Sect. 3) associated with the choice of the transport model (NEEM EU and Vostok holes). The transport configuration case is defined in terms of convection layer depth z conv and eddy diffusivity D eddy threshold. Note that the presented RMSD values depend on the calibration of the data and should be considered from a relative perspective. Recent updates in the calibration process explain the slight differences with the results presented in Buizert et al. (2012). and minor isotope (see e.g. Landais et al., 2006). The accuracy of the model for these species is quantified by the difference between the location of the modelled gravitational slope and the location of the gravitational slope observed on δ 15 N 2 (observed at 30 m). Using these indicators, we test different transport configurations involving the convective layer depth as well as the effective diffusivity definitions. The different transport configurations, illustrated in Fig. 1 and Table 1, are denoted with the index case.
The convective zone is introduced by considering an absence of gravitational settling, which is expressed in the model with w α = 0. The underlying assumption is that, as in the atmosphere, turbulent transport is fast enough to prevent the occurrence of gravitational settling. As turbulent transport (related to e.g., surface wind and pressure variations) is not explicitly represented in firn models, the convective layer thickness z conv remains constant and is imposed in order to fit the behaviour of inert gases (i.e., 4 m at NEEM according to Buizert et al., 2012; and 13 m at Vostok according to Bender et al., 1994). This allows for simulating the (isothermal) behaviour of δ 15 N, δ 86 Kr and δ 40 Ar in firn (see Fig. 1 for NEEM; δ 15 N at Vostok is illustrated in Sect. 4). It has a limited impact on reference gases at NEEM, whereas it significantly affects the RMSD at Vostok where the convective zone is much deeper.
The first test (case 1) suppresses gravitational settling in the convective zone, but does not include an eddy diffusivity. The RMSD is less than one for both sites and a flattening effect is observed on the inert gases in the convective zone. However, the calculated diffusivity in near-surface firn exceeds the CO 2 molecular diffusion in free air by 50 % at NEEM and 120 % at Vostok. The fact that molecular diffusion can not exceed the one in free air is then introduced (cases 2 and 3), possibly corrected to take into account the effect of upper firn porosity (case 4, with the c f factor). The diffusivity is, thus, set to D eddy when it exceeds the specified threshold according to Eq. (15). When the threshold depth for gravitational settling (z conv ) is not modified (case 2), inert gases are not affected and the impact on reference gases (RMSD) is small. By contrast, modifying the convective layer depth to make it equal to the diffusivity threshold depth calculated by the model (z conv = z eddy = 10.6 m at NEEM, case 3) more significantly increases the RMSD and shifts the gravitational slope of inert gases too deeply in the firn. Finally, using c f reduces the offsets observed on isotopes in the diffusive region and deep firn (reference case, case 4, with the imposed convective layer depth for gravitation and with the eddy diffusivity region starting at the surface and ending at z eddy = 18.8 m at NEEM, 39.8 m at Vostok). At NEEM, the values of z conv (4 m) and z eddy are, thus, very different. Other models obtain similar results: Fig. 5 in Buizert et al. (2012) shows that in all models, the effective diffusivity exceeds c F D CO 2 ,air (1.07 × 10 −5 m 2 s −1 ) in the ∼ 15-20 m range. Comparing the diffusivity profiles for the different test cases (top-left subplot in Fig. 1), most of the differences occur in the upper firn (down to 20 m below the surface), where the diffusivity is largely increased.
The relative importance of the advective and diffusive transport on CO 2 molecular flux computed with Eq. (12) is presented on Fig. 2, where DE08 is added in the comparison as a site with a large accumulation rate, no convective zone and an important advective flux. Figure 2 depicts the effects of advection and diffusion (molecular and eddy) on the trace gas fluxes, calculated respectively with ρ o α f (v + w air ) and ρ o α f ( w α − w α ) (using D eddy = 0 to model molecular diffusion and D CO 2 (z) = D CO 2 ,air = 0 for eddy diffusion). The ratio of advective to diffusive flux is quantified with the Péclet number, which we compute (choosing a unitary characteristic length) as: As w α is a function of ρ o α , ρ o α is first computed using the trace gas atmospheric scenario as input to the forward model (the date is chosen as 1990 to compute consistent fluxes for the different sites). On Fig. 2, the CO 2 fluxes illustrate the large variability of the transport processes depending on the site considered: mostly advective (below 30 m depth) at DE08, diffusive with a dominant advective contribution in the deep firn at NEEM, and diffusive down to the close-off depth at Vostok. The contribution of the eddy flux in the diffusion term is negligible at DE08 (only in the first metre), larger than molecular diffusion in the top 6 m at NEEM and has a strong impact at Vostok (larger than molecular diffusion down to 30 m). The Péclet number shows the transition between advection and diffusion (which occurs where Pe = 1). It quickly increases in the deep firn where the reduced diffusivity and the occlusion of air in bubbles play a major role on trace gas transport. The depth profile of Pe is illustrated for our eleven sites on Fig. 8 and is further discussed in Sect. 4.1. Several issues associated with discretisation are detailed in the Supplement, along with the definition of the numerical schemes. The main conclusions of this analysis are that Lax-Wendroff and first-order upwind methods can both be used safely for the discretisation of first order spacederivatives, and successive approximations of the solution with z = {0.8, 0.4, 0.2} m do not induce numerical instability, the scheme stabilising between z = 0.4 and 0.2 m (0.2 m is, thus, retained for the reference step size, a larger value inducing numerical errors). Implicit time discretisation with a sampling time of t s = 1 week provides a good tradeoff between numerical accuracy and simulation time (similar results are obtained with t s = 1 day while those obtained with t s = 1 month differ slightly). An explicit scheme with z = 0.2 m requires a computation time almost 500 times larger than the implicit model and is, thus, not suited for the optimisation procedures described in the next section, which involve multiple calls of the forward model. A mixed implicit-explicit scheme (Crank-Nicholson) provides similar results as the implicit scheme, but induced numerical instabilities in specific tests with the inverse scenario model.

Preliminary conclusions and perspectives
The forward model sequence of computations can be summarised as follows. Given the firn density profile, the full close-off depth and the site average temperature, we calculate the ice density as well as the total and open porosities from Eqs. (2), (3) and (6), respectively. The firn sinking velocity is then obtained from Eq. (1a) and the air velocity from Eq. (9). The trace gas concentration is computed by solving: in the convective layer (from the surface to z conv ) and Eq. (14) below z conv (diffusive and bubble closure regions). The boundary conditions are indicated in Eq. (14). Calculating the trace gas concentrations requires a diffusivity profile (including both CO 2 and eddy diffusivities; their computation is discussed in the next section), CO 2 diffusion coefficient in free air and the relative diffusion coefficient of the trace gas with respect to CO 2 in free air. The effective diffusivity is obtained from Eq. (15). The eddy diffusivity depth is automatically computed as the depth at which the effective diffusivity is r α c f D CO 2 ,air and is deeper than the convective layer (defined as the depth range where gravitational fractionation does not occur).
Modifying the advective-diffusive model allowed us to introduce the effect of gravity (non-isobaric condition) in the classical molecular diffusion framework. We could introduce this phenomenon despite the lack of information on permeability by considering the diffusive flow induced by a deviation from the steady-state (hydrostatic) equilibrium and modelled with Fick's law. While the resulting diffusive flux would be equivalent to the one proposed by Schwander (1989) if all the other model components were the same, our analysis leads to a better understanding of the underlying assumptions. Detailing the steps and assumptions used to introduce gravitational settling is also of interest in comparison with poromechanics, in which the gravity effect is generally included with Darcy's law (e.g., see Coussy, 2003). We further propose a new way of calculating D eddy in the upper firn based on the fact that molecular diffusivity in firn cannot exceed c f D CO 2 ,air .
The transport of gases with constant atmospheric concentration is well represented by the model, both in the diffusive region and in the firn-ice transition. Our D eddy approach of the convective zone, based on setting a maximum value for molecular diffusivity leads to fairly similar results as in other models (Buizert et al., 2012, , Fig. 5). A depth threshold for gravitational fractionation (z conv ) still needs to be imposed a priori. An alternative could be to introduce a variable surface pressure (statistically representative of atmospheric variations). The pressure variations due to wind stress on Atmos. Chem. Phys., 12, 11465-11483 a non-flat surface are more difficult to represent, but may strongly affect sites undergoing a high surface rugosity such as in areas of megadunes (Severinghaus et al., 2010). The steady-state approach to describe the non-isobaric property of firn can be used in other fields involving trace gas transport in porous media (with applications in soil sciences, porous catalysts or gas nuclear reactors, as mentioned by Webb and Pruess, 2003). It has the advantage of improving the advective-diffusive model without requiring information on permeability and to involve simple expressions for the diffusivity term.

Inverse diffusivity model
Inverse algorithms calculating firn diffusivity using firn-air observations of trace gases with well-known histories by minimisation techniques have been proposed in previous works by Rommelaere et al. (1997), Trudinger et al. (2002) and more recently by Buizert et al. (2012) in an intercomparison perspective. The forward model described in the previous section is robust with respect to diffusivity variations and is linear in terms of the trace gas concentration in open pores. It can, thus, be used in inverse diffusivity calculations. The nonlinear effects on the optimisation are due to the multiplication of the diffusivity by the trace gas concentration in the diffusive flux term of Eq. (14). Such effect is referred to as a bilinearity and has a limited impact on the problem of multiple solutions if the linearisation error is small (which is expected to be the case from the preliminary results for CH 4 at Dome C described in Witrant and Martinerie, 2010). Concerning the input data, the atmospheric scenarios and firn air measurements have to be consistent in terms of calibration scale and uncertainty estimates (as discussed in Martinerie et al., 2009;Buizert et al., 2012).

Formulation of the optimisation problem
The objective of the inverse diffusivity model is to determine the diffusivity profile that minimises the difference between the measured and the modelled mixing ratios. The diffusivity of each gas is parameterised in terms of the eddy and CO 2 diffusivities according to Eq. (15). The firn diffusive property is expressed in terms of the effective diffusivity D eff = {D eddy + c f D CO 2 ,air ; D CO 2 }. The optimisation (or inverse) problem is formulated for multiple gases as the distributed nonlinear least squares problem: where σ j,i is the standard deviation of the measured mixing ratio m j,i of trace gas j at the measurement location i, RMSD j is the root-mean-square deviation and N g is the number of trace gases considered. The notation "arg min" refers to the optimal diffusivity profile D * eff that minimises the cost function described by the sum term, D eff being the argument (free variable) of the optimisation problem. The final time t F corresponds to the measurement date. The time dependency is omitted for air since only hydrostatic equilibrium is considered and the gas concentrations ρ o α are obtained by solving the formard model up to t F . Additional physical properties of the firn diffusivity (D > 0 and ∂D/∂z < 0, i.e., diffusivity is monotonically decreasing with depth, which results in a negative gradient) are used as extra constraints to improve the convergence of the optimal solution. An analysis of this optimisation problem and solution was proposed by Witrant and Martinerie (2010) involving a linearisation of the diffusive flux in a partial differential equations framework.
The resulting inverse diffusivity model is presented in Fig. 3. First, for a given atmospheric scenario and D eff , the forward model provides the trace gas concentration distribution at final time ρ o j (z i , t F , D eff ). The firn air measurements m j,i and associated standard deviations are then introduced to compute the RMSD for each gas according to Eq. (19). The optimal diffusivity profile is finally obtained by iteratively updating D eff to minimise the resulting global error according to Eq. (18).

Optimisation algorithm and uncertainties
The minimisation problem is solved in a vector-valued approach. The vector containing the squared errors at the measurement depths for each gas is provided to the algorithm. The error vector is used to compute a preconditioned conjugate gradient (computed numerically using small variations on the diffusivity profile at each depth). The subspace trust-region method based on the interior-reflective Newton method (trust-region-reflective algorithm) described by Li (1994, 1996) then determines the modified diffusivity profile for the next iteration.
The convergence rate is improved by starting from a coarse grid (e.g., z = 0.8 m) and increasing the resolution once a minimum is reached. This method also decreases the impact of local minima and reduces computation time. The non-uniqueness of the solution (unavoidable considering our limited dataset) is illustrated in the results for each borehole presented in the Supplement by setting two different initial diffusivity profiles: one at zero and another set with a parameterised function. The differences in the final diffusivities are mainly limited to the upper firn and have no significant impact on the gas mixing ratios.
Specific care has to be taken estimating experimental uncertainties (σ i,j ), as they impact directly on the optimal solution. Two main sources of uncertainties are considered: -uncertainties on atmospheric trend scenarios σ scen , expressed as depth-dependent profiles by simulating their impact on concentrations in firn using maximum/minimum scenarios as input to the forward model (more details on this approach are provided in Buizert et al., 2012).
Note that we neglected the uncertainties of the measurement depths; justified by the fact that the gas is collected from a wider volume. For NEEM we used the standard deviation computation proposed by Buizert et al. (2012) while we applied a similar, but slightly different method (see Supplement) for other sites.

Reference and evaluation gases
We distinguish between reference gases (used to constrain the inverse diffusivity model) and evaluation gases (which allow us to discuss specific aspects of the physical model).
Among the evaluation gases, the isotopic composition of gases with constant atmospheric mixing ratios (such as δ 40 Ar, δ 86 Kr or δ 15 N) are of specific interest as they are not affected by uncertainties on the atmospheric scenarios. Therefore, they allow us to evaluate the steady-state behaviour of the forward model. However, due to their relatively low sensitivity to the diffusivity profile (gravitational fractionation is the main observed phenomenon in the diffusive region), they provide only weak constraints on the diffusivity profiles. This was observed at NEEM where adding δ 15 N as a reference gas did not significantly change the diffusivity profile or decrease the RMSD (with respect to our simulations with the nine other reference gases). This observation is in line with the system identification theory (see, e.g. Ljung, 1999), which states that the system inputs (atmospheric scenarios and firn air measurements in our case) have to be "sufficiently rich" (in terms of frequency content) to excitate all the model modes for a proper parameter identification. In contrast to our other reference gases where the frequency content in atmospheric scenarios partially compensates the limited depth resolution of firn data, the data associated with inert gases may not be informative enough. Therefore, we chose not to use inert gases in our diffusivity optimisation.
The relative weight of each gas used for the reference NEEM simulation is illustrated by the measurement signal to noise ratios (m j,i /σ j,i ) in Fig. 4. Although it does not account for the fact that gases which have steep concentration changes in the firn bring more information (e.g. Trudinger et al., 2012) the signal-to-noise ratio is closely related to the minimised cost function and illustrates the cost of the misfit of the firn data points. For NEEM EU, CH 4 is dominant in the upper 10 m while CO 2 has a major effect in the diffusive region. The reduced weight of CO 2 in the upper firn is due to the under-sampling of its seasonal variations . In the lock-in region 14 CO 2 contributes significantly (along with CO 2 and CH 4 ). Compared to other firn air models presented by , our (LGGE-GIPSA) model fits methane better than other species such as CH 3 CCl 3 and 14 CO 2 in their peak regions.
The final agreement between the firn model and specific data points is evaluated with the contribution of each measurement in the RMSD. The results depicted on the right panel of Fig. 4 show the overall coherency of the model (only two points with a significantly higher cost) and the increased difficulty to obtain a coherent matching in the convective and lock-in regions.

Sensitivity to the use of multiple gases at NEEM
In order to characterise the efficiency and limitations of the multi-gas approach, two tests are performed on the NEEM EU site. The first test, presented in Fig. 5, shows the effect of diffusivities optimised based on single gas measurements used on other gases (each reference gas is used in turn to compute the diffusivity). While a very good match is obtained for the gas corresponding to the computed diffusivity, a large divergence in model solutions appears both on the diffusivity and on the mixing ratios of other gases. This observation confirms that the inverse diffusivity problem is underconstrained for any single gas. It is interesting to note that gases with flat diffusivity profiles (e.g.,  in the upper layers of the firn) do not provide useful constraints on diffusivities in the associated region, as their mixing ratios are particularly invariant.
In the second test (Fig. 6), one of the reference gases is removed from the optimisation in turn (the diffusivity, thus, results from 8 gases out of 9). The robustness of the result (its ability to depict accurately the behaviour of several gases) is strongly improved, with very limited variations of the reconstructed diffusivity in the diffusive region of the firn. The Fig. 4. Signal-to-noise ratio of firn air measurements (left) and data weights in the minimised cost function (right) at NEEM (EU hole) for CO 2 (blue), CH 4 (green), SF 6 (red), HFC-134a (orange), CH 3 CCl 3 (violet), 14 CO 2 (turquoise), CFC-11 (pink), CFC-12 (yellow) and . The cost function distribution versus depth is presented for the reference case (one colour per gas).
remaining inconsistencies can possibly be induced by the 1-D hypothesis, the accuracy of gas diffusion coefficients, the accuracy of data uncertainties and the unmodelled transport phenomena. The upper region (top 15 m) still appears to be underconstrained as relatively large variations of the diffusivity do not induce significant variations of the gas mixing ratios. It is interesting to note that CH 4 and 14 CO 2 have an enhanced impact on the deep firn diffusivity. As these gases have a highly transient behaviour in the lock-in region, they bring more information to the multigas model than the other gases for this region.
Similar conclusions can be drawn on the NEEM US borehole (results provided in the Supplement), where we can also observe that using two gases already brings a significant improvement in comparison with a diffusivity profile constrained by a single gas.

Impulse response and the inverse scenario objective at NEEM
Age distributions constitute a model output of first interest for atmospheric scenario reconstruction, as discussed by Rommelaere et al. (1997) and Trudinger et al. (2002). The sensitivity of age distributions to the firn diffusivity, thus, provides an evaluation criterion for atmospheric trend reconstructions. Age distributions obtained from the impulse response of the model are presented for CO 2 at NEEM (EU hole) in Fig. 7. A convective zone test case is provided by comparing a case where the eddy diffusivity is set to zero (case 1) with the reference case (case 4), discussed in Sect. 2.4. While this clearly has a significant impact on the upper firn, the difference is relatively small below 60 m, especially near the age of maximum probability (maximum gain amplitude). Another test case is provided by setting the diffusivity to zero in the deep firn (when it becomes smaller than 0.05 m 2 yr −1 , corresponding to the lock-in region). The surface response is not affected, but the responses below 70 m are significantly modified (13 % difference in the peak gain at 78.6 m). As most of the information for past history reconstruction older than 20 yr is contained below 65 m, diffusivity values in the upper firn are less important than in the diffusive region and in the lock-in region.

Trace gas transport at Arctic and Antarctic sites
In addition to NEEM-EU and NEEM-US, eleven firn air pumping operations previously modelled (Rommelaere et al., 1997;Fabre et al., 2000;Sowers et al., 2005;Faïn et al., 2009;Martinerie et al., 2009) with the single-gas diffusivity minimisation algorithm of Rommelaere et al. (1997) were simulated using our new multi-gas method. In this section, we try to decipher links between site specific behaviour and major climatic characteristics: temperature and snow accumulation rate. Gaining knowledge about the relationship between firn physics and climate is important for the interpretation of trace gas records in ice cores (see e.g. Landais et al., 2006). Firn models with a tuned diffusivity profile can be directly used only for present-day firn air pumping operations for which depth profiles of trace gases in firn are available. Modelling past firn profiles requires to approximate their diffusivities with a simple parameterisation dependent on climatic parameters, which is provided at the end of this section.

δ 15 N and gravitational behaviour
In this section, our forward model results using the optimal diffusivity profile are compared with δ 15 N data at the studied drill sites (Fig. 8) in order to evaluate the capability of the model to predict the location of the lock-in zone. Thus, we focus mainly on the deepest part of the δ 15 N profiles. In the upper firn, δ 15 N is affected by thermal diffusion driven by the seasonal variations of temperature (Severinghaus et al., 2001), a process not represented in our model.
At intermediate depths, model results and δ 15 N data (within uncertainties) primarily follow the gravitational slope. However, the model deviates from the gravitational slope at the highest accumulation rate sites (DE08 and Devon Island). This confirms the observation of Trudinger et al. (1997) that faster firn sinking (advection) prevents δ 15 N from reaching gravitational equilibrium. It also implies that the advective velocity associated with firn sinking plays a major role in our capability to capture the deviation from the gravitational slope. Battle et al. (1996) introduced the name of lock-in zone, defined as the deep firn region where δ 15 N becomes constant instead of being enriched by gravitation. They associated the absence of gravitational enrichment to an absence of diffusive mixing. However, the effective firn diffusivity in the upper part of this region is still significant (see Sect. 4.2), as observed by Buizert et al. (2012) for the NEEM firn. The LIZ is defined as the depth range between the end of δ 15 N gravitational fractionation (referred to as the lock-in depth LID) and the full bubble close-off depth (z F , f = 0). Within Atmos. Chem. Phys., 12, 11465-11483 Multiple gases inverse diffusivity model for NEEM (EU hole): the diffusivities (reduced to CO 2 ) are computed with 8 reference gases out of 9 (one gas removed for each test). The removed gases are CO 2 (blue), CH 4 (green), SF 6 (red), HFC-134a (orange), CH 3 CCl 3 (violet), 14 CO 2 (turquoise), CFC-11 (pink), CFC-12 (yellow), CFC-113 (brown) and none (black, reference).
experimental uncertainties, the model captures well the width of the LIZ. Another important observation on Fig. 8 is that the width of the LIZ increases when snow accumulation increases (the main physical characteristics of the 13 modelled drill holes are summarised in Table 2). The most arid Antarctic plateau sites (Dome C and Vostok) virtually have no δ 15 N plateau (or δ 15 N defined LIZ). The same behaviour occurs at Megadunes (Severinghaus et al., 2010). Thus, the location and width of a LIZ derived from δ 15 N does not represent the bubble close-off zone at low-accumulation sites (otherwise the bub-ble close-off zone would have to be very narrow and occur at depths where air can hardly be pumped out of the firn). As glacial periods are more arid than interglacial periods, the decreasing width of the δ 15 N LIZ associated with the decrease of snow accumulation rate can generate a bias when trying to infer past bubble close-off depths from δ 15 N data in ice cores. In Table 2, we compare three definitions of the lock-in depth: the LID derived from δ 15 N (depth at which the model results flatten), the approximate depths at which a slope break is observed in the CO 2 and/or CH 4 depth profiles (called LID gas , see e.g., Fig. 1) and the depth at which Table 2. Physical characteristics of modelled firn air pumping sites. Devon Island, Summit, North Grip and NEEM are Arctic sites; other sites are Antarctic sites (DML stands for Dronning Maud Land). Major climatic parameters: snow accumulation rate (in centimetres water equivalent per year) and temperature are indicated in columns three and four. Columns 5 to 7 provide quality indicators for the diffusivity profiles: N gas is the number of reference gases used for constraining diffusivity, the RMSD column provides the root-mean-square deviation between model results and data for the reference gases (the number between parentheses is the difference between the maximum and minimum RMSD calculated with N gas − 1). Columns 8 to 14 provide indicators of the top and bottom heights of the LIZ: the last sampling (or measurement) depth, the model's full COD (open porosity is zero, corresponding to z F in Eq. 4), LID gas is the approximate depth at which a slope break is observed in the CO 2 and/or CH 4 depth profile (see e.g., Fig. 1), LID D is the depth at which D CO 2 ≤ 1 m 2 yr −1 , δ 15 N LID is the depth at which δ 15 N fractionation stops, LID Pe is the depth at which the Péclet number becomes larger than 5, 10 % COD and 50 % COD are the depths at which the closed/total porosity ratios are 10 % and 50 %, respectively. The last column provides the mean age at 50 % COD.

Site
Drill Accu. Temp. N gas RMSD Last meas z F LID gas LID D δ 15 N LID LID Pe 10 % COD 50 % COD 50 % COD year (cm w eq)  Fig. 7. Contribution of air of different ages at given depths (impulse response of the forward model for CO 2 at NEEM EU hole) using the optimal reference diffusivity (case 4, "-"), without eddy diffusivity (case 1, "---"), and with eddy diffusivity, but cancelling the diffusivity when D < 0.05 m 2 yr −1 ("· · · "). the closed/total porosity ratio is 10 % (called 10 % COD), using the modified parameterisation from Goujon et al. (2003) described in Sect. 2.1. δ 15 N LID and LID gas are consistent within uncertainties (2-3 m) except at Dome C and Vostok where they differ by 5-6 m. Figure 9 further illustrates the relationship between the lock-in zone width and the accu-mulation rate. It shows estimated ranges of the LIZ widths defined using δ 15 N (in black) and a reference gas (in grey). The increase of the LIZ width with accumulation rate mostly occurs at accumulations lower than 10 centimetres, and δ 15 N LIZ widths are frequently lower than gas LIZ widths. The difference between δ 15 N and gas LIZ widths is highest at the very arid sites (Dome C and Vostok). By contrast, the 10 % COD is deeper than δ 15 N LID by 2 to 6 m at sites with 10 cm water equivalent per year or more snow accumulation, 10 % COD and δ 15 N LID are approximately equal at sites with accumulation rates of about 7 cm water eq. yr −1 , and 10 % COD is about 10 metres higher than δ 15 N at the lowest accumulation sites (< 4 cm). If the widths of the LIZ defined by δ 15 N or LID gas appear related to the snow accumulation rate, the width defined by the 10 % COD shows a more erratic behaviour (it varies between approximately 6 and 11 m without a clear link to the accumulation rate). The LID can also be related to the transition from diffusive to advective transport in deep firn by calculating LID Pe , the depth at which the Péclet number defined in Eq. (16) reaches an upper bound. The effect of a high value of tortuosity preventing gravitational enrichment in deep firn is quantified with LID D , the depth where the effective diffusivity becomes smaller than a lower bound. Choosing the bounds Pe = 5 for LID Pe and D CO 2 = 1 m 2 yr −1 for LID D , we obtain the results presented in Fig. 10 and Table 2. LID D provides an estimate that is very close to LID gas , which reflects the model capability to capture observed changes in the slope of the reference gases into the diffusivity profile. LID Pe provides an estimate that is closer to the δ 15 N LID at Atmos. Chem. Phys., 12, 11465-11483 Nitrogen (δ 15 N of N 2 , left scale) as a function of depth at different sites: measured ("*") versus modelled (blue "-") isotopic composition and comparison with the gravity slope ("---"), along with the Péclet number calculated from the concentration profile of CO 2 (green, right scale). δ 15 N stands for δ 15 N of N 2 (in ‰). The firn air sampling sites are shown by order of decreasing snow accumulation rate (see Table 2). Experimental data for NEEM are corrected for the effect of thermal diffusion and are average results for the NEEM EU and NEEM US drill holes . low accumulation sites. LID Pe is less accurate at predicting δ 15 N LID than LID D at predicting LID gas , it shows important deviations at Siple, Berkner and DE08. On the other hand, the δ 15 N defined LIZ is clearly different from the bubble close-off zone, at least for low accumulation rate sites. The different LID estimates are very consistent at the best-constrained site (NEEM). As the model's open to closed porosity parameterisation is scaled to the firn density profile and based on data at only three sites (Goujon et al., 2003), a better physical characterisation of the deep firn (lock-in and close-off zones) in future firn air sampling programmes would help improving the ice-core related understanding of firn physics.

Diffusivity profiles
Table 2 also provides quality indicators related to the diffusivity profile. The primary diagnostic is the root-mean-square deviation between model results and data for all gases used. Fig. 9. Ranges of δ 15 N defined (in black) and reference gas defined (in grey) lock-in zone widths versus accumulation rate. The plotted ranges take into account the large uncertainty on the full close-off depth (z F ). Uncertainties on z F were estimated using available information about depths at which air could or could not be extracted from the firn and/or firn air sampling resolution.
However, it is more difficult to simultaneously fit a large number of species than a small number of species, thus, for example, the RMSD is higher for NEEM-EU (nine gases) than NEEM-US (three gases). A useful complementary indicator is the difference between the minimum and maximum RMSD calculated with N-1 gases, this difference is much smaller at NEEM-EU (0.06, the smallest of all sites) than for NEEM-US (0.18). Note that all of our calculated mean RMSD values are lower than one (see Table 2). However, these low values are difficult to interpret in an absolute sense , as they can be due to either pessimistic estimates of the data uncertainties or an overfitting of the noise in the data. An in depth discussion of this point can be found in Trudinger et al. (2012). Figure 11 compares the diffusivity profiles obtained for our 13 target drill holes. The fit of reference gas data at each site is illustrated in the Supplement, together with sitespecific tests. In the upper firn (open porosity > 0.3) Vostok, North GRIP and NEEM-EU show higher diffusivities than other sites, in relation with the occurrence of specific convective zone processes. Such processes are also indicated by δ 15 N for Vostok (Bender et al., 1994) (13 m deep convective zone), NEEM-EU  (∼4 m deep convective zone), and North GRIP where nearly constant mixing ratios of several gases are observed in the top ∼8 m (see Supplement). On the other hand, Devon Island shows reduced diffusivities in most of the firn due to the presence of about 150 refrozen melt layers (see e.g. Martinerie et al., 2009). Although these heterogeneities are not modelled, the multigas diffusivity constraint brought a strong improvement at this site (see Supplement). Diffusivities at other sites than Devon Island are remarkably consistent for open porosities in the range 0.1-0.3. Understanding the deep firn results requires assessing where diffusivity becomes negligible. Sensitivity tests were performed where firn diffusivity was set to zero when reaching different threshold values. Setting to zero the values of the firn diffusivity scaled by the CO 2 diffusion coefficient in free air (D eff /D air ) which are lower than 10 −3 has a significant impact on reference gas mixing ratios in all the LIZ, and on the δ 15 N LID, whereas setting D eff /D air values lower than 10 −4 to zero has a very small influence on the results. Thus, the large range of model values below D eff /D air = 10 −4 are insignificant. For values of D eff /D air higher than about 10 −3 (corresponding to D eff > 0.5 m 2 yr −1 , as D air ≈ 500 m 2 yr −1 ), the decrease of diffusivity with decreasing open/closed porosity ratios is slowest for the lowest accumulation rate sites (Vostok and Dome C), shows intermediate values at South Pole and DML and is fast at other sites (the somewhat slower decrease at DE08 may be related to uncertainties due to the low sampling resolution in the LIZ for this early drill site and/or an insufficiently constrained diffusivity due to the dominant effect of advection). The higher deep firn diffusivities at low accumulation rate sites are consistent with faster transport and, thus, the accumulation rate dependent behaviour of δ 15 N illustrated in the previous section. We, furthermore, evaluated gas ages at the depth where f/ = 0.5 (50 % COD age in Table 2). Very young deep firn ages calculated at Dome C and Vostok also support a fast gas transport in deep firn at these sites. Note that very young gas ages were also calculated with the Rommelaere et al. (1997) model for Dome C (Martinerie et al., 2009). For other sites, the behaviour of the 50 % COD age is more complex. It shows a general trend of increasing ages when accumulation decreases. Some exceptions may be related to the variable thickness of the firns (e.g., younger ages at DML than South Pole). Other exceptions have no obvious explanation, such as the distinctly younger ages at Summit than other Arctic sites (North GRIP and NEEM), while these three sites have similar temperatures and accumulation rates and Summit is the deepest firn.
A scaling law has been estimated based on the above diffusivities at the different sites, including the site temperature T and atmospheric pressure P atm . We obtain the approximation: where T 0 and P 0 are the reference temperature (273 K) and pressure (101 325 Pa). The LIZ is included by defining a threshold value D thr on D eff , with the corresponding depth z thr . Below z thr , D eff is approximated with a sigmoid curve as: D eff (z) ≈ D thr − 10 −2 1 + e 50(z−(z F −z thr )/2)/z F + 10 −2 A reasonable approximation of the LID is obtained by distinguishing the sites with low accumulation rates from the sites with high accumulation rates as: D thr = 1 m 2 yr −1 if a accu > 10 cm w eq., D thr = 100 m 2 yr −1 otherwise. A comparison between this estimation and the optimum diffusivity at each site is provided in the Supplement, along with the resulting CO 2 and CH 4 mixing ratio estimates.

Conclusions
The formulation of trace gas transport in firn (forward model) was revised in a poromechanics framework leading to a quasi steady-state approach for including the effect of gravity. The effect of diffusion, specifically the turbulent behaviour of gases (eddy flows) in the upper firn, was investigated in relation with the convective layer behaviour.
We introduced a new multi-gas constrained inverse model for diffusivity calculation, based on nonlinear least squares minimisation. The inverse problem formulation, optimisation algorithm and choice of the reference gases were discussed. The available measurements do not allow for the inverse problem to be fully constrained. However, the multi-gas approach greatly improves the efficiency of the inverse diffusivity model. Evaluation tools such as the signal-to-noise ratio of gas measurements, space distribution of the modelling errors and age distributions provided insights on the inverse model efficiency and limitations. Despite two strong constraints in the diffusivity calculation: diffusivity decreases with depth and no eddy term is considered in the lock-in zone, this new model showed very good performance for NEEM in a recent inter-comparison study . The proposed new model was applied to thirteen firn air pumping sites. The analysis of stable isotope ratios of molecular nitrogen (δ 15 N) for the different sites allowed us to validate the capability of our model to derive the lock-in depth location from mixing ratio measurements of reference gases. This depth threshold could be associated with a reduced diffusivity value (high tortuosity) and a large Péclet number (dominant advective flow), and qualitatively related to the accumulation rate. The reduced width of δ 15 N flat part (which defines the lock-in zone) at low accumulation rate sites also appears associated with higher diffusivities in the upper bubble close-off region (10 −2 > D eff /D air > 10 −3 ). For the most arid sites, where no flat part is observed, anomalously young mean gas ages are also obtained. Taking into account the accumulation rate dependency of the width of δ 15 N-defined lock-in zone may lead to a better understanding of the difference between ice age and gas age in ice cores. A scaling law of the diffusivity profile including a LID-specific behaviour was proposed to model trace gas transport in past firns.
Future work investigating the deep-firn physics (e.g., tomography) and providing a temporal tracking of the upper firn behaviour would be highly valuable to improve the understanding of these least well constrained parts of the firn. Evaluating the firn permeability depth-profile would also be of major interest to capture transient effects induced by seasonality and fast atmospheric pressure variations.