How well do different tracers constrain the ﬁrn diffusivity proﬁle? Supplementary material

Abstract. Firn air transport models are used to interpret measurements of the composition of air in firn and bubbles trapped in ice in order to reconstruct past atmospheric composition. The diffusivity profile in the firn is usually calibrated by comparing modelled and measured concentrations for tracers with known atmospheric history. However, in most cases this is an under-determined inverse problem, often with multiple solutions giving an adequate fit to the data (this is known as equifinality). Here we describe a method to estimate the firn diffusivity profile that allows multiple solutions to be identified, in order to quantify the uncertainty in diffusivity due to equifinality. We then look at how well different combinations of tracers constrain the firn diffusivity profile. Tracers with rapid atmospheric variations like CH 3 CCl 3 , HFCs and 14 CO 2 are most useful for constraining molecular diffusivity, while &delta: 15 N 2 is useful for constraining parameters related to convective mixing near the surface. When errors in the observations are small and Gaussian, three carefully selected tracers are able to constrain the molecular diffusivity profile well with minimal equifinality. However, with realistic data errors or additional processes to constrain, there is benefit to including as many tracers as possible to reduce the uncertainties. We calculate CO 2 age distributions and their spectral widths with uncertainties for five firn sites (NEEM, DE08-2, DSSW20K, South Pole 1995 and South Pole 2001) with quite different characteristics and tracers available for calibration. We recommend moving away from the use of a firn model with one calibrated parameter set to infer atmospheric histories, and instead suggest using multiple parameter sets, preferably with multiple representations of uncertain processes, to assist in quantification of the uncertainties.


Ice conservation and coordinate systems
Unit volume has proportion s of air and 1 − s of ice, since mass (per unit volume) is (1 − s) × ρ ice = ρ, giving total porosity s(z) = 1 − ρ(z)/ρ ice (1) The mass of ice (per unit area) above z is and vertical velocity of ice (relative to the surface) is We wish to use a moving coordinate system, so care is needed in addressing transformations from a fixed coordinate system. This is done in terms of a generic function F (.,.), with the need to consider that as coordinates change, the functional relation will also change, and that the notation should reflect this.
Here we define a number of coordinate systems, along with the transformations and notation for functions. We are interested in the final coordinate system (y,t ), and the others are given only to show how (y,t ) relates to (z,t).
The τ coordinate can be linearly rescaled to give units in terms of distance by multiplying by any fixed velocity. The most obvious choice is surface velocity v(0) = A/ρ(0) ξ,t . Where ξ = τ − t − a with a an arbitrary constant. This coordinate moves with the ice. We use t = t, with the distinction between t and t helping clarify the meaning of the various partial derivatives.
y,t . This is a re-scaled version of ξ, scaled by v(0) = A/ρ(0). giving and Note that rather than using d dt and ∂ ∂t notation to distinguish derivatives in moving and fixed coordinates, our t and t do this. We will express some equations in terms of (t , z) -this allows easier comparison between equations in fixed and moving coordinates, and is more convenient for implementation.
If we use a 'tilde' notation to indicate functions defined in the (y,t ) coordinates then the defining relation is: F (z,t) =F (y(z,t),t (t)) For the present analysis, we assume that the ice properties are constant with time (the case of a melt layer moving with the ice is an exception that will be considered later), so ice properties are functions of z only. For any quantity G that is independent of t, the relations above give: and In the y,t coordinates, z dependence can be mapped onto either rates of change with respect to t or y.

Pore space
The total porosity s is partitioned into open porosity f and closed porosity b, as: Moving with the ice, closed porosity changes due to (i) new trapping (which increases b with depth) and (ii) bubble compression (which decreases b with depth). Open porosity changes due to (i) new trapping (which decreases f with depth) and (ii) compression of channels (which decreases f with depth). By writing b(z) = r(z)s(z) and f (z) = (1 − r(z))s(z), we can separate the effects of trapping and compression on f and b: ∂b ∂z = r ∂s ∂z + s ∂r ∂z (14) Like Rommelaere et al. (1997) and Severinghaus and Battle (2006), we assume that the volumetric compression acts equally on open and trapped pore space so that we interpret the term involving ∂s ∂z as compression and s ∂r ∂z as trapping. Figure S1 shows open, closed and total porosity at NEEM, Northern Greenland, with density and open porosity values from the NEEM intercomparison (Buizert et al., 2012). It also shows their derivatives with depth, and these split according to Eqs. (14) and (15), indicating where compression and trapping have most influence on f and b. For example, closed porosity is altered most by trapping, but there is also an effect due to compression shown by the long dashes in

Air conservation in pore space
Like Trudinger et al. (1997) and Rommelaere et al. (1997), we assume that air in open channels, C air , can be described by the barometric equation Fig. S1. a) Open (f ), closed (b) and total (s) porosity at NEEM, b) fraction of total porosity that is closed r, c) ∂s ∂z , ∂f ∂z and ∂b ∂z , d) ∂r ∂z , e) ∂f ∂z split into the contributions from compression (1 − r) ∂s ∂z and trapping −s ∂r ∂z , f) ∂b ∂z split into the contributions from compression r ∂s ∂z and trapping s ∂r ∂z , g) vertical velocities (in m/yr) of ice (v), air (w) and upward flux of air due to compression (−u), h) air pressure in the trapped bubbles (solid) and open firn (dotted).
In fixed coordinates: and Moving with the ice, we defineC air (y(z,t),t (t)) so that, as above Air in open channels in the ice is compressed as the ice moves down to regions of lower total porosity, and this leads to an upward flux of air (relative to the ice) to maintain the balance described by the barometric equation. We define this flux, denoted by u, as positive downwards, so u is negative.
Air is also trapped into new bubbles, and we assume that a new bubble formed at depth z has the same pressure (when formed) as air in open pores at that depth.
In fixed coordinates, air is carried down with the ice at velocity v(z) but is also expelled upwards to give net velocity w(z) = u(z) + v(z). Air conservation in open channels in fixed coordinates can be written as (Rommelaere et al., 1997) The left side of the equation is zero, the first term on the right describes flux divergence and the second term describes trapping into new bubbles. Using the product rule on the first term on the right hand side, with Eq. (20) for ∂C air ∂z and dividing by C air and f we get This equation can be solved for w with the boundary condition w = v where f goes to zero. We then calculate u from u = w − v. We found this easier than solving the equivalent moving coordinate equation for u, particularly at the boundary where f goes to zero, and the equations are equivalent. Rommelaere et al. (1997) pointed out that in fixed ice coordinates, the net vertical flux of air downwards, w(z), balances trapping into bubbles. Figure S1g shows the vertical velocities v, u and w calculated for NEEM.

Air conservation in bubbles
We need to keep track of the air pressure in the trapped bubbles to determine the change in mole fraction of tracers as new bubbles are formed. After the first bubbles are formed, they move down with the ice and compression causes the closed porosity to decrease and the air pressure to increase. Then as more bubbles are trapped with the same air pressure as in the open pores at that depth, the closed porosity (of bubbles in a layer of ice) increases and the air pressure decreases. Air conservation in trapped bubbles can be determined from (Rommelaere et al., 1997) which we will need later, and We can solve this equation using finite difference representation of the derivatives with the boundary condition at the surface C b air (0) = C air (0) Equation (27) can also be derived by considering the compression and trapping alternately for each time step, as in Rommelaere et al. (1997). Figure S1h shows the calculated air pressure in bubbles and open firn at NEEM. The calculated air pressure in the bubbles depends on the form of the closed porosity variation with depth. At NEEM we have used the Goujon et al. (2003) parameterisation from Buizert et al. (2012) where closed porosity increases gradually from the surface but is negligible until around 60 m. A spline fit to closed porosity measurements such as that used in Trudinger et al. (1997) gives quite a different variation with depth for the calculated C b air , but the value when trapping stops is also around 1.4 × C air (0).

Trace gas conservation in open pore space
The concentration of a generic trace gas, C x (z,t), (in mol m −3 ) is modelled (in fixed coordinates) in terms of the flux through open firn channels, J x (z,t), given by The first term describes molecular diffusion, the second term gives the settling due to gravity. The third term describes an 'eddy' diffusion, as introduced by Severinghaus et al. (2001), which is the same for all trace gases and primarily includes mixing near the surface due to convection. Note that the eddy diffusion term is parameterised as the deviation from a hydrostatic gradient (J. Severinghaus, pers. comm. 2011). The last term is the total downward flux of air due to advection of ice w = v + u.
The conservation equation (in fixed coordinates) is defined in terms of the full volume where the gas flow per unit area is f (z)J(z,t) and the amount of gas per unit volume is The term involving λ is radioactive decay and the last term describes bubble trapping.
We can derive the equivalent equation in moving coordinates using w = u + v and Eq. (10) which can be rearranged to give where the second term on the left side of Eq. (32) has cancelled with a term on the right side, and the last term on the right side of Eq. (33) is the other term from the product rule applied to the left side of Eq. (32). This is the equation for trace gas concentration in moving coordinates. This equation differs from the previous version of the CSIRO firn model  in that the flux now includes eddy diffusivity and the upward flux of air due to compression u, and the mass balance equation includes three extra terms. These equations give results as concentrations, C x (z,t), in mol m −3 . However, as mole fraction in dry air is more commonly measured than concentration, we wish to write the diffusion and mass conservation equations in terms of mole fraction, c x (z,t), where C air (z) is described by the barometric equation (Eq. (18)), giving and in fixed coordinates. giving and from we can substitute in Eqs. (35) and (36), then use Eqs. (9) and (10), as follows We replace C x by c x C air in Eq. (33), expand and a number of terms will cancel, to leave Note that a number of terms cancelled using Eq. (24). If we define the flux of tracer x in terms of mole fraction in moving coordinates as then we have These equations now only differ from those in Trudinger et al. (1997) by the addition of the eddy diffusion term, and the upward compression flux −u ∂cx ∂z . The other differences that we saw in the concentration equations disappeared when we converted to mole fraction. Notice that the bubble trapping term has disappeared, as it doesn't directly affect the mole fraction in the open firn.
The factor v(0) v in front of the spatial derivatives appears because we are expressing the derivatives in terms of y rather than z. Instead of solving for the derivatives in terms of y, we can solve these equations but evaluate the derivatives in terms of the physical depth z, using

Trace gases in trapped bubbles
We model the trace gas mole fraction for both open and closed pores in all model layers. As soon as closed porosity increases above zero, the model will start to calculate the mole fraction in the trapped bubbles, although if this occurs gradually in the shallow firn (such as in the Goujon et al. (2003) parameterisation), the modelled trapped mole fraction can be associated with a negligible amount of air. Once diffusion has stopped, air is locked into the channels in a particular piece of ice, and the mole fraction of any trace gas in that model layer stays constant (apart from radioactive decay) as the air is progressively trapped into bubbles. We model the mole fraction in the closed pore space, c b x , by considering both the compression of closed pore space and trapping of air into new bubbles. Compression increases the air pressure and trace gas concentration in the bubbles, but not the trace gas mole fraction. As new bubbles are formed we combine the previous trace gas mole fraction in the closed pores with new bubbles.
Conservation of mass gives leading to It we use the quotient rule on c b This weights the trace gas mole fraction already in old bubbles and that trapped into new bubbles by air content, taking into account compression of the old bubbles that has already occurred. The old CSIRO model weighted by porosity, and didn't take into account compression. The difference between the two methods is actually very small, but taking into account compression is expected to be a more accurate representation of reality. We have assumed that there is no fractionation due to bubble trapping, which may occur for smaller molecules than those considered here (Huber et al., 2006;Battle et al., 2011).

Implementation
To derive an equation (in moving coordinates) that we can solve with the implicit time stepping, we can expand Eqs. (45) and (46) and collect terms involving ∂ 2 cx ∂z 2 , ∂cx ∂z and c x to get As already mentioned, the variable y was useful to ensure that the equations were derived correctly in moving coordinates, but it is more convenient to evaluate the derivatives in terms of z. However, it is very important to account for the fact that layers have different thicknesses in terms of z. When evaluating derivatives in space, we use and where c k , c k−1 and c k+1 are the concentration in the layers k, k − 1 and k + 1, and h 1 = z k − z k−1 , h 2 = z k+1 − z k and h 3 = (h 1 + h 2 )/2 are distances between layer centres. If we want to remove the fractionating effect of gravity (for comparison with measurements that have been corrected for gravity), we can set M x = M air . Isotopic ratios are modelled in the firn by modelling each of the isotopomers as separate tracers with their own atmospheric history as a mole fraction, then combining the model output of the tracers at each depth to give the isotopic ratio.
Model layers correspond to equal mass of ice per unit area (in kg m −2 ), denoted W and usually a fraction or multiple of the annual accumulation. Each new layer accumulates in time τ = W/A. The time interval τ is divided into an integer number of timesteps. Initially the surface layer has mass W , and each timestep mass is added to this layer and the layers below it are moved downwards, until the surface layer reaches mass 2W . At this point the coordinate system is relabeled, such that the values of quantities in layer k are put into layer k + 1 and the surface layer is reset to mass W . At each model time step we need to recalculate the depth of each layer and the corresponding density, open porosity and diffusivity before solving the model equations.
The model is usually run with a time step dt=0.01 yr. In order to speed up the GA calculation, each time a new parameter set is tested the firn model is initially run with dt=0.5 yr until a few years before the end of the calculation when we change to 0.01 yr to the end (to capture the variation causes by the seasonal cycle). If Φ for this parameter set is less than 0.4 (for NEEM) above the upper bound of solutions we are saving, the model is rerun with dt=0.01 yr. The difference in Φ between these two cases (at NEEM) varies between about zero and 0.3. The initial fast run is about 20 times faster than the rerun. If a wide prior parameter range is used, there can be many parameter sets tested that are a very bad fit to the data, so running first with a large time step offers a significant time saving compared to running all cases with dt=0.01 yr.

Convective mixing
We have two options for modelling convective mixing near the surface. The first is a well-mixed layer, similar to that described by Trudinger et al. (2002). In the new version of the model, the model layers start at the depth of the well-mixed layer, as in Rommelaere et al. (1997), where this depth can be tuned by the GA along with diffusivity. The second option for modelling convective mixing involves using an exponentially decreasing eddy diffusion following Severinghaus et al. (2001), where the two parameters describing eddy diffusion (magnitude and length scale) can be tuned by the GA.

Melt Layers
The ice structure at DE08-2 shows a melt layer at 8.7 m below the surface. Trudinger et al. (1997) found that the agreement between modelled and measured tracers at DE08-2 (SF 6 in particular) was significantly improved by including a melt layer that originated at the surface in the 1989-90 summer and moved with the ice with a reduction of diffusive flux of about 80%. Sofen (2007) also considered a melt layer at Summit. A melt layer moving with the ice is a departure from our assumption that the ice properties at a particular depth are constant with time. For simplicity, we assume that the melt layer affects only molecular diffusion but not any of the other physical properties of the ice (open or closed porosity, bubble trapping, convective mixing/eddy diffusion, air pressure or air flow). When air was collected at DE08-2 in 1993, the melt layer was at 8.7m, so was too shallow to have affected bubble trapping yet. Thus, the assumption that the melt layer has affected only molecular diffusion seems appropriate for DE08-2 but may not be suitable for other sites, particularly for a deeper melt layer. We include the melt layer in the model for DE08-2 by replacing the diffusivity D X (z ML ) at the layer boundary corresponding to ice that fell as snow in 1989.77 (i.e. depth z ML ) by µD X (z ML ) where 0 ≤ µ ≤ 1 and µ is the degree to which the melt layer has reduced molecular diffusion. Although the model equations were derived assuming that the ice properties were constant with time, there were no equations involving time derivatives of diffusivity, so adding time variation to D X does not pose a problem. Time variation of open or closed porosity would need much more consideration. Although multiplying D X (z ML ) by µ is equivalent to multiplying the flux J X (z ML ) by the same value, as was done in Trudinger et al. (1997), there are other differences between the old and new versions of the model (upward flux of air due to compression in the new version and flux smoothing in the old version) that mean that the value of µ here may give a different effect compared to the equivalent parameter in the old model. Here we find that a reduction of diffusion of 89% gives the optimal match to DE08-2 observations.

Atmospheric histories
We use the atmospheric concentration histories for the high latitude northern hemisphere from Buizert et al. (2012), shown in Fig. S2. These were compiled from direct atmospheric measurements, firn/ice core measurements from Law Dome in Antarctica, tree-ring data and emissions-based estimates. Atmospheric records were at monthly resolution starting in the year 1800.
The high latitude southern hemisphere atmospheric histories for CO 2 and CH 4 are from the Law Dome ice core records (Etheridge et al., 1996(Etheridge et al., , 1998MacFarling Meure et al., 2006) and Cape Grim direct atmospheric measurements (Francey et al. (2010) and references therein; Rigby et al. (2008) and references therein). For SF 6 , CFCs, HFC-134a and methyl chloroform we use histories based on AGAGE measurements and emission-based model results from Martinerie et al. (2009), consistent with the northern hemisphere records. For 14 CO 2 the southern hemisphere history is based on measurements from tree rings and from Wellington, New Zealand (Manning and Melhuish, 1994), and HCFC-141b is based on measurements at Cape Grim (Cape Grim Air Archive (unpublished data, P. Krummel per-sonal communication) followed by direct atmospheric measurements (O'Doherty et al., 2004)). Table S1. DE08-2 firn measurements (Etheridge et al., 1996;Trudinger et al., 1997;Levchenko et al., 1997;Etheridge et al., 1998). CO2 and CH4 are given in the WMOX2007 and NOAA04 scales, respectively. SF6 is given in the University of Heidelberg SF6 scale. Measurements of ∆ 14 CO2 in permil have been converted to 14 CO2 using CO2 and δ 13 CO2 measured on these firn samples. We use the same uncertainties for all measurements of each tracer, with values given on the last line of the

DE08-2 and DSSW20K firn measurements
The firn measurements that we use for DE08-2 (Etheridge et al., 1996;Trudinger et al., 1997;Levchenko et al., 1997;Etheridge et al., 1998) are given in Table S1, and for DSSW20K (Smith et al., 2000;Sturrock et al., 2002;Trudinger et al., 2002) in Table S2. For DSSW20K, we chose not to use measurements of SF 6 , CFCs and methyl chloroform below 45 m in our model calibration, because their atmospheric histories (for the southern hemisphere) prior to 1978 are based on emissions estimates rather than atmospheric measurements. These measurements are given in the table in brackets. The uncertainties used are shown on the last row of the tables. HCFC-141b uses uncertainties of 0.1 ppt for the upper two observations and 0.05 ppt for the rest.

NEEM Synthetic B calculations
The 'Synthetic B' observations have the same true concentration as Synthetic A, but observations correspond exactly to the types and measurement depths for the EU borehole  in the NEEM intercomparison (between 15 and 23 measurements of each tracer). Three types of error are added to the synthetic observations: (1) random error with zero mean and sd = 0.5 % of the range, (2) systematic offset added to each tracer, where the magnitude of the offset (a constant value for each tracer) is a random number with sd = 1 % of the range for each tracer, and (3) systematic error that increases linearly with depth from zero at the surface to a value at 80 m that is a random number with sd = 2 % of the range for each tracer. Depending on the random values generated, some tracers will have larger systematic errors added than others. The data uncertainty used to calculate Φ B is σ i = R× (0.005 2 + 0.01 2 + (0.02 · z i /80) 2 ) where R is the range for each tracer described above. This case is more like reality, where systematic errors probably dominate, than Synthetic A. The uncertainties on the NEEM observations from Buizert et al. (2012) were often around 1-2 % of the range for each tracer, with some larger values for various reasons, so this case is similar to the real NEEM case in terms of the general error characteristics (but not in terms of the specific errors on each tracer). Figure S3 shows the diffusivity profiles and CH 3 CCl 3 determined by calibrating with these observations. With a lowest value of Φ B of 0.8 for Subset Ten, we consider solutions with Φ B up to 0.84 corresponding to the 68% confidence level for all Synthetic B subsets. We also show results for Table S2. DSSW20K firn measurements (Smith et al., 2000;Sturrock et al., 2002;Trudinger et al., 2002). CO2 and CH4 are given in the WMOX2007 and NOAA04 scales, respectively. SF6 is given in the University of Heidelberg SF6 scale. Halocarbon measurements from Sturrock et al. (2002) are given on the SIO2005 scale. Measurements of ∆ 14 CO2 in permil have been converted to 14 CO2 using CO2 and δ 13 CO2 measured on these firn samples. The lower measurements of SF6, CFC-11, CFC-12, CFC-113 and CH3CCl3 in brackets are not used for calibration, because their atmospheric histories are based on emission estimates rather than atmospheric measurements. The same uncertainties are used for all measurements of each tracer, with values given on the last line of the  Figure S4 shows the synthetic diagnostic scenarios used in Buizert et al. (2012) to compare different aspects of the firn models. As well as results of the other models from Buizert et al. (2012), we show our results for the new CSIRO model with the NEEMTenEddy case for a range up to Φ N = 0.92 (the range covered by the models in Buizert et al. (2012)). The scenarios are described in detail in Buizert et al. (2012), and only very briefly here. The first scenario (Fig. S4a) compares diffusive fractionation for a hypothetical monotonic CO 2 increase. The range from the different models in Buizert et al. (2012), with molecular and/or eddy diffusion in the lock-in zone, is larger than our range that includes cases with only molecular diffusion. The second scenario (Fig. S4b) shows attenuation of a 15yr period sinusoidal CO 2 with depth. Here our range is larger than the range from Buizert et al. (2012). The third scenario (Fig. S4c) shows gravitational enrichment for gas X with a very small relative diffusion coefficient γ X = 0.025. Our range is smaller than the range from different models. The fourth scenario, mean age of gas Y with advective transport only (γ Y = 0, Fig. S4d

Inferring relative diffusion coefficients
As well as there being uncertainty in the CO 2 diffusivity profile, there is also uncertainty in the relative diffusion coefficients, γ X = D X /D CO2 . Previous firn studies have used quite a wide range of relative diffusion coefficients. For example, some values of γ CH4 used in past firn modelling studies include 1.35 for Summit (Schwander et al., 1993), 1.29 for DE08-2 , 1.415 for DE08-2 (Rommelaere et al., 1997; Martinerie et al., 2009) and the value used here of 1.367 (Buizert et al., 2012). For γ SF6 there has been 0.582 , 0.621 (Martinerie et al., 2009) and 0.554 used here from Buizert et al. (2012). This variation corresponds to roughly ±5% variation around the middle of the range; we expect less than 1% variation due to differences in temperature. The different estimates of γ X were based on measurements, data compilations and empirical equations from various sources (Andrussow et al., 1969;Marrero and Mason, 1972;Fuller et al., 1966;Chen and Othmer, 1962;Lugg, 1968). The values of γ X used in Buizert et al. (2012) were based on a consistent set of diffusivity measurements from Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 with estimated uncertainty of about ±2%. Due to the different estimates available for relative diffusion coefficients, Trudinger et al. (2002) tuned γ SF6 and γ HFC134a at DSSW20K, obtaining values 0.628 and 0.614, respectively. Trudinger (2000) found best agreement to South  Pole CH 4 with γ CH4 of 1.42. Butler et al. (1999) also adjusted γ CFC-11 by 10 % to fit observations better. Here we investigate tuning the relative diffusion coefficients with the GA using both the pseudo and real observations. We repeated the Synthetic A and B Subset Ten calculations estimating the relative diffusion coefficients, γ X , of seven tracers in addition to the diffusivity profile and well-mixed layer depth. We allowed a range for the diffusion coefficients of ±10 % of the true value. Sensitivity tests showed that δ 15 N 2 varies very little when γ15 N2 is varied within this range, so we did not try to estimate it; we also did not try to estimate γ14 CO2 , but included observations of these two tracers and used their 'true' values of γ X . Figures S5 and S6 show scatter plots of diffusion coefficients against Φ for all solutions tested by the GA that had Φ A < 1.25 for Synthetic A and Φ B < 1.0 for Synthetic B. The clustering of points into horizontal lines is due to the way the GA algorithm works, re-taining solutions with low Φ and mutating or breeding them. Synthetic A has a steeper rise in Φ for a similar change in relative diffusion coefficient (i.e. a narrower minimum) than Synthetic B. Synthetic A also has more sparse coverage of solutions than Synthetic B in these plots -many more solutions are generated by the GA in Synthetic B in our range of interest than in Synthetic A, despite the same specifications for the GA. This was a feature of scatter plots of all parameters in the synthetic calculations -for Synthetic A they were sparse but Synthetic B they were dense.
For most tracers the Synthetic A case gives a greater mismatch from the noisy concentration observations and from the true concentrations than Subset Ten with the true relative diffusion coefficients. This is probably because we are now estimating a larger number of parameters (23), so it has become harder for the GA to locate the true solution. The solution with the lowest Φ A has estimated relative diffusion co- Table S3. Weighted RMS mismatch of modelled concentrations from observations and from the truth for each tracer for each of the five Subsets using Synthetic B observations. Modelled concentrations are for the case with lowest ΦB for each experiment. Numbers in bold are for tracers that were fitted for that particular Subset. As well as RMS mismatch for each tracer, we also show the RMS mismatch for all ten tracers together (All), and the RMS mismatch for only those tracers that were fitted in that experiment (Fitted). Data uncertainties are used for weights in the cost function, and we use Synthetic A data uncertainties as weights for the mismatch from the truth for Synthetic B concentrations to allow better comparison with Synthetic A results. The first column of numbers shows the RMS mismatch of the true solution from the noisy observations weighted by the data uncertainties.

Mismatch from observations, ΦB
Mismatch from truth, ΦAt Truth Two Three Four Five Ten Two Three Four Five Ten    efficients that are between 0.9 % too low and 2.4 % too high. The scatter plots show how well resolved each estimate is. For Synthetic B, the best solution for most tracers is closer to the noisy observations but further from the truth than the case with known diffusion coefficients. The best case has diffusion coefficients between 2.1 % too low and 8.4 % too high. CFC-11 in particular has an estimated diffusion coefficient that is very high, and it is the tracer with the largest systematic errors added (apart from δ 15 N 2 and 14 CO 2 for which we're not estimating γ X ), see the first column in Table S3.
The Synthetic A calculation does quite a good job at estimating γ X , particularly as we are now estimating 23 parameters at once. Most of the best estimates of γ X are slightly high, but there is only one tracer tying them together (CO 2 ), so we might expect a slightly lower CO 2 diffusivity versus depth profile compared to the truth to compensate. The Synthetic B results are more concerning, as the estimates are further from the true γ X and often outside the ±2 % uncertainty range quoted for the Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 values. We repeated the calibrations for NEEM, DE08-2 and DSSW20K with relative diffusion coefficients estimated along with the other parameters. For NEEM, we recalibrated Subsets Ten and TenEddy estimating relative diffusion coefficients. The best cases had Φ N significantly lower than the values obtained with fixed diffusion coefficients (Φ N = 0.66 compared to 0.74 with fixed γ X for TenEddy). There was most reduction in mismatch for SF 6 and CFC-12; there appears to be a slight calibration bias between NEEM CFC-12 in the upper firn and the atmospheric history, and this has probably affected the result. The best case had relative diffusion coefficients ranging from 11 % higher to 8 % lower than the values given in Buizert et al. (2012). CH 4 and CFC-12 had the greatest change, with the best Φ N obtained with lower values of relative diffusion coefficient, as well as a higher value for SF 6 . Figure S7 shows the results for all three sites, adjusted to apply to a temperature of 244.25 K. We started with a range of ±10 % of the Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 values (shown by the error bars in Fig. S7), but if the GA preferred a solution at the boundary we extended the range and reran the calculation. To the left of each of the error bars, we show various published estimates (measured and empirical) of γ X for eight tracers. To the right of the error bars are the values we determined using the GA for the five firn sites (with NEEM included twice, using each of the methods for convective mixing near the surface). These values give a better fit to the observations than using the Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002) γ X , however we do not know whether the different values of γ X are compensating for other errors, such as in the atmospheric history or missing or incorrectly modelled processes in the firn model. Scatter plots of relative diffusion coefficients as a function of Φ in each case show a clear preference for particular values of relative diffusion coefficient, and often (although not always) a fairly steep increase in Φ as you move away from that value.
In some cases the values estimated by the GA are close to the Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 values, but in other cases they are significantly different. It is difficult to trust these results after the Synthetic B calculations, if we assume that real observations could suffer from similar systematic error. If Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 are correct that their diffusion coefficients are accurate to ±2 %, then their values are more reliable than our calibrated estimates, certainly at present.  S7. Relative diffusion coefficients, γX , for eight tracers. The error bars indicate ±10 % around the Matsunaga et al. (1993) values (or Chen and Othmer (1962 for CH3CCl3 and Fuller et al. (1966) for HCFC-141b where Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 values are not available). The letters to the left of the error bars are measured or empirical estimates of γX , with M = Matsunaga et al. (1993), F = Fuller et al. (1966, C = Chen and Othmer (1962), L = Lugg (1968) and m = Marrero and Mason (1972). The letters to the right are estimates from our GA calibrations for different firn sites: N = NEEM with a well mixed layer, n = NEEM with exponential eddy diffusion, D = DE08-2 and W = DSSW20K. All values are adjusted to apply to temperature of 244.25 K. The right axis applies to CH4 and the left axis applies to all other tracers.
However, an uncertainty of even ±2 % for the diffusion coefficients estimated by Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002, and probably higher for other tracers not estimated by Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 such as CH 3 CCl 3 and HCFC-141b, will contribute to the overall uncertainty in the inferred diffusivity profile and therefore the uncertainty in atmospheric reconstructions. One way to account for this uncertainty would be to include relative diffusion coefficients in the GA calibration, as we have already done, but with a range of ±2 % (or greater for some tracers) for the purpose of taking account of their uncertainty rather than to improve the estimates of their values. We would probably prefer to use as our best solution a case with the Matsunaga et al. (1993Matsunaga et al. ( , 1998Matsunaga et al. ( , 2002 values of γ X , even though other values of γ X might give a lower Φ, but we could easily include some additional members in our representative ensemble that have these other values of γ X to reflect this component of the uncertainty. As firn measurements and models become more accurate this may become a relatively more important contribution to the overall uncertainty.

Dispersion in the lock-in zone
Dispersion in the lock-in zone has recently been included in a number of firn models to improve the fit to observations (Severinghaus et al., 2010;Buizert et al., 2012). The dispersive mixing term is mathematically identical to the eddy diffusion term, with a transport flux that is equal for all tracers and their isotopologues. Dispersion therefore gives no gravitational separation with depth. Severinghaus et al. (2010) and Buizert et al. (2012) discuss the physical reasons behind the idea of dispersive mixing in the lock-in zone, including upward mixing of old air due to the increase of pressure in summer layers before bubbles are formed. In Buizert et al. (2012), all six firn models required nonzero diffusivity of the order of 0.1 m 2 y −1 in the lock-in zone, and this was parameterised in different ways using either molecular diffusion, eddy diffusion (dispersion) or a mixture of both. Four of the models had dispersive mixing in the lock-in zone to improve the fit to the slow-diffusing gases such as CFC-113 (Severinghaus, 2012), or because it reduced the peak in 14 CO 2 closer to observations while leaving δ 15 N 2 constant below the lock-in depth. The other two models, including the CSIRO model, used only molecular diffusion in the lock-in zone. There was no obvious difference in Φ between models using the different types of diffusion in the lock-in zone. Our new results for NEEM show a very slight increase in δ 15 N 2 with depth down to about 77 m (about 0.004 ‰) that is about as consistent with the measured δ 15 N 2 as constant levels.
The different treatment of mixing in the lock-in zone in Buizert et al. (2012) was believed to have caused significant variation in estimates of diffusive fractionation affecting trace gas isotopic ratios, such as δ 13 CO 2 and δ 13 CH 4 . More work is definitely needed to understand the processes occurring in the lock-in zone and to confirm whether dispersion occurs. However, because dispersion may be important for diffusive fractionation in particular, here we show some preliminary calculations that add dispersion to our GA calibration to see how well it can be constrained and how much the spread of solutions increases. These calculations are a start to accounting for model error in our uncertainty estimates.
As a first test we added dispersive mixing to calibrations with NEEM synthetic observations, estimating dispersion for two cases: one with dispersion in the forward run to generate the synthetic data, and one without (both cases had errors like Synthetic A). We specified dispersion as a function of open porosity, f , prescribed as a non-symmetric cosine function (i.e. a cosine function between −π and π where the peak could be shifted left or right by scaling each side differently in the horizontal direction). We describe this function with four parameters, as shown in Fig. S8: the maximum height of the function D max , the porosity value at the left (−π) side of the function f min , the width of the function (in terms of porosity), f wid , and the location of the maximum of the function as a fraction of the width, f frac (0.5 would put the peak in the middle, giving a symmetric function). For the open porosity range f min < f < f min + f frac f wid , the equation is and for the range This form is not as general as our function for molecular diffusion, but serves as a starting point to allow some dispersion in the lock-in zone with a few parameters that can be tuned. For dispersion in the forward run we used the cosine function with D max = 0.1 m 2 yr −1 , f min = 0.002, f wid = 0.17 and f frac = 0.65. In both cases we included the 4 parameters describing the dispersion peak in the GA calibration as well as the usual molecular diffusivity and the depth of the well-mixed layer, to see whether we could resolve dispersion in either case.
We found a difference in the results for the two cases, as shown in Fig. S9. The results on the left show the case with dispersion in the forward run used to generate the pseudo observations. Mismatch is denoted Φ Ad , and we consider solutions with Φ Ad < 1.25. Results on the right are without dispersion in the forward run, and we consider solutions with Φ A < 1.25. The best solution in the case with dispersion is similar to the true dispersion function used, and the difference of the best estimate from the true dispersion is offset by a similar but opposite difference between the true molecular diffusivity and the best estimate. The true molecular diffusion is unfortunately outside the envelope of solutions with Φ Ad < 1.25 for part of the range. In the case without dispersion, there are some accepted solutions with dispersion above 60 m, where molecular diffusion dominates, and there are some narrow dispersion peaks accepted below about 70 m, but the best estimate of molecular diffusion is very close to the truth.
We were therefore able to recover roughly the correct dispersion profile in both cases using observations with small, Gaussian noise. While there was information in these observations to distinguish between molecular and dispersive mixing in the lock-in zone, observations with realistic errors would be expected to reduce the ability of the calibration algorithm to distinguish between molecular diffusion and dispersive transport. Note that the case with dispersion had both molecular diffusivity and dispersion in the lock-in zone, and we could distinguish them to some extent. They might be easier to distinguish if there was one or the other.
With 11 tracers for calibration including 14 CO 2 , DSSW20K is another site where we might be able to resolve dispersion in the lock-in zone. We are also interested to see whether allowing dispersion in the lock-in zone means that the calibration no longer prefers exponentiallydecreasing eddy diffusion that extends down to the lock-in zone. We ran two additional calibrations using DSSW20K observations and dispersion in the lock-in zone, one with the well mixed layer and the other with exponentially-decreasing eddy diffusion for convective mixing near the surface (not cut off at 30 m). We allowed the magnitude of the dispersion to vary up to 10 m 2 yr −1 .
The lowest values of Φ W for these two additional DSSW20K cases were 0.95 with the well-mixed layer for convection (previously 1.12 without dispersion) and 0.91 with eddy diffusion (previously 0.92). The well-mixed layer cases are not quite as good as the eddy diffusion cases. Cases allowing dispersion in the lock-in zone fit observations more closely than cases without dispersion -this is not just a case of fitting observations more closely when there are more parameters to tune, as our original case with exponential eddy diffusion extending into the lock-in zone had no extra parameters. It really seems like dispersion in the lock-in zone allows a better match to DSSW20K observations. However, as none of the cases fit the peak in 14 CO 2 very well, this leaves us a little cautious of this conclusion, and we would like further confirmation along with better understanding of the actual processes involved and estimates of the likely magnitude of dispersion. The GA preferred dispersion close to 10 m 2 yr −1 for the case with the well mixed layer and around 2 m 2 yr −1 for the other case (which already had eddy diffusion from the exponential function). These levels are considerably higher than the values used in Buizert et al. (2012) of around 0.1 m 2 yr −1 and Severinghaus et al. (2010) of order 0.01 m 2 yr −1 . When we allow dispersion in the lock-in zone, the exponentially-decreasing eddy diffusion still extends into the lock-in zone, but to a much lesser degree.
We selected representative subsets of 20 parameter sets in each of the five DSSW20K cases, and calculated the isotopic diffusion correction. This correction is used to correct observations of δ 13 CO 2 for fractionation in the firn due to the different rates of diffusion of the isotopes . For this calculation, we used spline fits to CO 2 and δ 13 CO 2 measurements from Law Dome and Cape Grim. Figure S10 shows the calculated diffusion correction with depth for cases without any eddy diffusion in the lock-in zone in the upper panel, and all cases in the lower panel. All solutions have Φ W < 1.17, which corresponds to the 68% con- Fig. S10. a) Diffusion correction for δ 13 CO2 in the DSSW20K firn. The black line shows our preferred case for calibration with exponentially-decreasing eddy diffusion cut off at 30 m and no dispersion in the lock-in zone. The blue lines show results for 19 representative solutions for this case. Red lines show the diffusion correction for 20 solutions with the well-mixed layer. b) Same as a, but adding 20 solutions with exponentially-decreasing eddy diffusion extending into the lock-in zone (purple), 20 solutions with exponentially-decreasing eddy diffusion cut off at 30 m and dispersion in the lock-in zone (orange), and 20 solutions with a wellmixed layer for convective mixing and dispersion in the lock-in zone (green). All solutions have ΦW < 1.17. fidence level for our preferred case. Including dispersion in the calculation has roughly doubled the diffusion correction range at its maximum.
Our calculations for DSSW20K show that unless it is specifically avoided by selection of parameter prior ranges, it is possible for the exponentially-decreasing eddy diffusion for convection near the surface to extend through the whole firn and have an influence well into the lock-in zone. The values of the parameters preferred by the GA for DSSW20K were probably chosen to suit both convection near the surface and dispersion in the lock-in zone, rather than just convection near the surface which was the intention. Our specification of molecular diffusion was chosen to avoid being prescriptive about the variation with depth, but the specification for eddy diffusion (both for convection near the surface and dispersion in the lock-in zone) is quite prescriptive. In the SIO model in Buizert et al. (2012), the balance between molecular diffusion and dispersion was specified with a single coefficient varying between 0 and 1, rather than trying to estimate a depth-varying dispersion profile that is poorly constrained.
Adding in the possibility of dispersion or eddy diffusion leads to greater equifinality than with molecular diffusion alone, but if evidence for its existence continues to grow it should be accounted for as part of the range of possibilities. There is much less chance of resolving dispersion at sites with fewer tracers for calibration. Tracers with γ X further from 1.0 or isotopic ratios are more likely to distinguish between molecular diffusion and dispersion than tracers with γ X near 1.0.