Numerical and analytical solutions for σ

In Iso-CFM, we follow the original time-stepping scheme of the CFM (Fig. 1). For each iteration j, using the output of CFM for dρ/dt, ρ and T we calculate the quantity dσ²/dt as:

(22)$${{\rm d}\sigma^2\over {\rm d}t} = 2 \left(D\lpar {t} \rpar - {\sigma^2\over \rho} \, {{\rm d}\rho\over {\rm d}t}\right)$$

In Figure 2, we plot the three terms of Eqn (22) as well as the diffusion length value that is the result of the numerical integration of Eqn (22). The forcing we use is T = 242 K and $A = 0.131 \, {\rm m\, a^{-1}}$ ice eq. and we implement the HLD model. As seen in Figure 2, the diffusion length signal is a result of two opposing processes, on one hand the positive contribution of the isotope firn diffusivity (2D(t), dot line in Fig. 2) and on the other hand the negative contribution of the densification process (term ${-2\sigma ^2\over \rho }{{\rm d}\rho \over {\rm d}t}$, dot-dash line in Fig. 2). The diffusivity term 2D(t) is positive at all times and decreases from the surface and until the close-off depth where it is equal to zero. Note how at depths below ≈ 30 m in our calculation the densification term dominates over the diffusion term and thus dσ²/dt is negative and the value of the diffusion length is decreasing. Note also the discontinuities on both rate terms at ≈ 15 m. These are due to the transition zone of the HLD model at the critical density of 550 kg m⁻³ described by the two activation energies in Eqns (26) and (27).

Fig. 2. Contribution of the diffusion and densification terms (Eqn 22) on the evolution of the diffusion length signal as a function of depth.

For the initial iteration of the model, we calculate the diffusion length profile using analytical equations derived from Eqn (8). Calculations based on these analytical equations have previously been shown in Johnsen and others (Reference Johnsen2000); Gkinis and others (Reference Gkinis, Simonsen, Buchardt, White and Vinther2014) and Holme and others (Reference Holme, Gkinis and Vinther2018). However, to our knowledge the equations themselves have not been published. Therefore, we present them here with a short description of their derivation. Rearrangement and substitution of variables in Eqn (8) results in:

(23)$${{\rm d}\sigma^2\over {\rm d}\rho} + {2\sigma^2\over \rho} = 2 {\left({{\rm d}\rho\over \rm{d}t}\right)}^{-1} D\lpar \rho\rpar .$$

The latter can be converted into the integral form:

(24)$$\sigma^2 \lpar \rho \rpar = {\, 1\, \over \rho^2}\int_{\rho_o}^{\rho}2\rho^2 {\left({{\rm d}\rho\over {\rm d}t}\right)}^{-1}\! D \lpar \rho \rpar \, {\rm d}\rho.$$

We use the densification rate parameterisation from Herron and Langway (Reference Herron and Langway1980):

(25)$${{\rm d}\rho\lpar z\rpar \over {\rm d}t} = k\lpar T\rpar A^{\vartheta} \lpar \rho_{{\rm ice}} - \rho\lpar z\rpar \rpar \comma\; $$

where k(T) is a temperature-dependent Arrhenius-type densification rate coefficient described by:

(26)$$k_0\lpar T\rpar = 0.011 \exp \left(- {10\, 160\over RT} \right)\comma\; \vartheta = 1\comma\; \; \; \; \rho\lt 550 \; {\rm kg\, m}^{-3}\comma\; $$

and

(27)$$\eqalign{k_1\lpar T\rpar = 0.575 \exp \left(- {21\, 400\over RT} \right)\comma\; \vartheta = 0.5\comma\; \; \; \; \rho\ge 550 \; {\rm kg\, m}^{-3}.}$$

Using the diffusivity coefficient from Johnsen and others (Reference Johnsen2000) (the quantities involved in the calculation of the diffusivity coefficient are described in detail in Section ‘Implementation of the firn diffusivity’) and substituting the fraction term in front of the parentheses with ${\rm \Xi}$

(28)$$D \lpar \rho \rpar = \overbrace{{mpD_{ai}\over RT\alpha_{i}}}^{\Xi} \left( 1 - {1.3\rho^2\over \rho_{\rm ice}^2} \right) \left( {1\over \rho} - {1\over \rho_{\rm ice}} \right) $$

we finally obtain the analytical equations for the diffusion length σ for the sections above and below the critical density ρ_c = 550 kg m⁻³:

(29)$$\sigma^2 \lpar \rho \lt \rho_{\rm c} \rpar = {{\rm \Xi}\over \rho^2 k_o A \rho_{\rm ice}} \left[\rho^2 - \rho_o^2 - {1.3\over 2\rho_{\rm ice}^2} \lpar \rho^4 - \rho_o^4 \rpar \right]$$

(30)$$\eqalign{\sigma^2 \lpar \rho \ge \rho_{\rm c} \rpar = {{\rm \Xi}\over \rho^2 k_1 A^{0.5} \rho_{\rm ice}} \left[\rho^2 - \rho_{\rm c}^2 - {1.3\over 2\rho_{\rm ice}^2} \lpar \rho^4 - \rho_{\rm c}^4 \rpar \right]\cr + {{\rm \Xi}\over \rho^2 k_o A \rho_{\rm ice}} \left[\rho_{\rm c}^2 - \rho_o^2 - {1.3\over 2\rho_{\rm ice}^2} \lpar \rho_{\rm c}^4 - \rho_o^4 \rpar \right]}$$

Iso-CFM calculates both the analytical solution of the diffusion length, as well as the time derivative dσ²/dt from Eqn (22).

Contour of analytical solutions for the diffusion length at the close-off density

Based on Eqn (30) we calculate the diffusion length value at the close-off density for ${\rm H}_2{}^{18}{\rm O}$ using a broad range of temperature-accumulation rate combinations. Hereafter we use the prime notation σ′₁₈ to denote the diffusion length value at the close-off depth. The calculation results in a contour plot (Fig. 3) that (a) allows the estimation of the expected diffusion length value for a number of ice coring sites and (b) provides an overview of the uncertainties expected in reconstructing temperatures based on a diffusion length value. It is apparent that an uncertainty of 1 cm for the diffusion length value σ′₁₈ results in different uncertainties in temperature. As seen in the contour plot, at the accumulation level of 0.1 $\hbox{m}\, \hbox{a}^{-1}$ ice eq. a σ′₁₈ uncertainty of 1 cm results in a temperature uncertainty of ≈ 5 K for temperatures around 220 K and of ≈ 2 K for temperatures around 240 K. This is an auxiliary plot to some of the diffusion length uncertainty calculations following in the paper.

Fig. 3. Contour map of analytical solutions for σ′₁₈ expressed in cm of firn eq. at the density of ρ _co. For all the calculations P = 0.7 Atm,  ρ_o = 350 kg m⁻³,  ρ _co = 804.3 kg m⁻³. The dash straight line represents the temperature–accumulation forcing of the steady-state experiments as described by Eqn (31).

Steady-state climate forcing

We run two experiments using steady-state climate forcing. For both experiments, the model is spun up for 1000 years, and the main model run lasts 400 years. The model domain is set to 300 m and heat diffusion is enabled. The first experiment (A), consists of 30 model runs, each with a unique temperature and accumulation rate pairing. Following previous studies (Johnsen and others, Reference Johnsen, Dahl Jensen, Dansgaard and Gundestrup1995; Dahl Jensen and others, Reference Dahl Jensen1998), we assume a logarithmic relationship between the forcing temperature and accumulation. The temperatures range linearly in the range 213–250 K, and the accumulation rates follow a natural logarithm relationship described as:

(31)$$\ln{A} = -21.492 + 0.0811\, T$$

Equation (31) describes a temperature–accumulation relationship that should be seen as qualitative only and not necessarily as representative of actual ice sheet conditions.

The temperature–accumulation forcing is shown in Figure 4. For experiment A we use annual time steps. Experiment B consists of seven runs for which we force the model using monthly time steps with a seasonal temperature cycle (Thompson, Reference Thompson1969; Lundin and others, Reference Lundin2017):

(32)$$T = T_{\rm{amp}}\lpar \cos\lpar 2\pi t\rpar + 0.3\cos\lpar 4\pi t\rpar \rpar $$

where T _amp = 10 K. Note that no seasonality is assumed for the accumulation, which remains constant at its mean annual value. The accumulation–rate temperature pairs for experiment B are shown as pink diamonds in Figure 4.

Fig. 4. Steady-state T and A forcing for experiments A and B. The pink markers indicate the mean forcing for experiment B for which temperature a seasonal cycle with an amplitude of 10 K is enabled.

For every run, diffusion length profiles for δ¹⁷O, δ¹⁸O and δD are calculated. In Figure 5 we present the diffusion length value for each isotopologue at the close-off depth (ρ _co = 804.3 kg m⁻³) expressed in m of firn (vertical distance) and symbolised with the primed notation as σ′_D, σ′₁₈ and σ′₁₇. We also include the calculations for the close-off depth. For the specific temperature–accumulation relationship we use here (Eqn 4), it is clearly seen that the effect of the increasing temperature is dominating by increasing the diffusion length value and shifting the close-off depth further up. For the first effect, the driving mechanism is the firn and air diffusivity as described in Eqns (9) and (13). The higher densification rates resulting in a shallower close-off depth can be explained by the temperature dependence of the densification rates which in the case of the HLD and HLS are given in Eqns (26) and (27). The BAR and GOU models use similar relationships although with different activation energies. The increase in accumulation rate results in a more effective ‘sealing’ between the layers thus reducing the isotopic gradients and hence hindering diffusion. Additionally, higher accumulation values result in a deeper close-off depth and as a result allow for a diffusion process that lasts longer in time. These two effects of the accumulation increase are not visible in this experiment due to the temperature dominating the process. For a better insight in the separate effects of temperature and accumulation on the value of the diffusion length, we refer the reader to Sections ‘Ramp experiments’, ‘Warming pulse experiments’ and Figures 12, 20 and 21.

Fig. 5. Results of model runs using steady-state forcing with and without temperature seasonal cycle enabled. Solid lines (markers) represent runs with the temperature seasonal cycle disabled (enabled). The prime notation for the diffusion length σ′ represents the value of the diffusion length at the close-off depth (ρ = 804.3 kg m⁻³) and expressed in m of firn.

Overall, the results of experiment A indicate a very good agreement between the HLD, HLS and BAR model whereas the GOU model deviates from the rest of the models in the range of higher forcing temperatures/accumulations. Specifically, as the temperature reaches values close to 250 K, the discrepancy is in the order of ≈ 1 cm for σ′₁₈ (similar for σ′₁₇ and σ′_D). As far as the seasonal temperature signal is concerned, we observe a minimal influence of the annual cycle on the results with an effect that is more profound for the case of the GOU model. In Figure 22 in the Appendix, we give a more detailed view of the temperature and diffusion length signals with and without seasonality.

The influence of the fractionation factor parameterisation

We investigate the influence of the parameterisation used for the fractionation factor $\alpha _{\rm s/v}^{i}$. by performing a test identical to experiment A and by using the various expressions presented here for $\alpha _{\rm s/v}^{i}$ (Eqns 14–18). As seen from Figure 6, the calculated discrepancies are minimal. A notable exception is that of the Ellehøj and others (Reference Ellehøj, Steen-Larsen, Johnsen and Madsen2013) parameterisation for deuterium that seems to yield a higher (although still not significant) deviation from the results obtained using the Merlivat and Nief (Reference Merlivat and Nief1967) and Lamb and others (Reference Lamb2017) relations. Based on the good agreement between the experiments of the two latter studies, we choose to use the Merlivat and Nief (Reference Merlivat and Nief1967) and Majoube (Reference Majoube1971) parameterisations throughout the paper.

Fig. 6. Estimation of the diffusion length value σ′ at the close-off depth (dash lines) and the firn diffusivity for the density of ρ _c = 550 kg m⁻³ (solid lines) with various parameterisations of the fractionation factor $\alpha _{\rm s/v}^{i}$. For the plots where there appears only one curve per parameter, the lines are visually indistinguishable from each other.

Diffusion length profiles

We present (Fig. 7) density and diffusion length profiles for the last time step of the model (year 400) for two sets of climate forcing regimes. The first one (type 1) is roughly representative of conditions of the East Antarctic Plateau (T ₁ = 222.66 K, $A_1 = 0.031 \, \hbox{m}\,\hbox{a}^{-1} \hbox{ice} \;\hbox{eq.}$) and the second (type 2) is more representative of conditions of ice coring sites on the Greenland ice sheet ($T_2 = 242 \, {\rm K}\comma\; \, A_2 = 0.131 \, \hbox {m\, a}^{-1} \hbox { ice\; eq.}$) (Table 1).

Fig. 7. Firn density and diffusion length profiles for type 1 (dash lines) and type 2 (solid lines) steady-state climate forcing. In (c) we also plot the annual layer thickness λ_A deduced from the density profile using the HLD model.

Table 1. Type 1 and type 2 climate forcing regimes

First, we note the good agreement between the H–L type models. The two dynamic implementations (HLD and HLS) and the analytical H–L model yield very similar results for all the diffusion length and density calculations. This is not a surprising result, considering that all three implementations are variations of the same model. However, as we show later, the dynamic response of these three implementations can be rather different. Second, there is a significant discrepancy between the GOU and the rest of the models for the type 2 regime (warmer conditions). This discrepancy is in line with the results of Section ‘Steady-state climate forcing’.

The lower diffusion length values of the GOU implementation for the type 2 forcing can be explained by the higher densification rates predicted by this model for the full span of the firn column. Similarly, the BAR model densifies faster than the HLS and HLD models for the upper part of the firn. However, for densities in the range ≈ 600–800 kg m⁻³, the H–L-type models predict higher densification rates. The net effect is that the three models yield very similar values for the diffusion length at the close-off depth. This indicates that both the upper ( < 550 kg m⁻³) and the lower stages ( > 550 kg m⁻³) of the densification play an important role in the diffusion process.

As seen in Figure 7, for type 1 conditions the diffusion length is higher than the annual layer thickness throughout the whole firn column, thus obliterating the annual cycle. For the higher accumulation rates of the type 2 climate forcing we should expect a better ‘sealing’ of the adjacent firn isotopic layers, effectively hindering diffusive transport. However, it is apparent that the temperature of the firn column plays a dominant role in the formation of the σ₁₈ signal resulting in higher diffusion length values.

Close-off density and the tortuosity factor

The value of the close-off density ρ_co signifies the point at which 1/τ → 0 and thus D = 0 (Eqn 9). At this point the diffusive fluxes seize and the process of vapour diffusion terminates. Commonly, in gas-diffusion firn studies the term ‘full close-off depth’ refers to the depth at which the open porosity of the firn reaches the value of zero and hence the bubbles are occluded (Martinerie and others, Reference Martinerie, Raynaud, Etheridge, Barnola and Mazaudier1992). Typical values for the full close-off density are ≈ 830 kg m⁻³. Despite the risk of confusing these two different terms we will use the term ‘close-off depth’ and ‘close-off density’ in order to be compatible with previous studies on water isotope firn diffusion.

The empirical form for 1/τ (Eqn 19) is based on laboratory measurements of the diffusivity. The majority of the measurements are performed for the CO₂ and O₂ gases using real firn core sections (Schwander and others, Reference Schwander, Stauffer and Sigg1988). In the study of Jean-Baptiste and others (Reference Jean-Baptiste, Jouzel, Stievenard and Ciais1998), diffusivity coefficients were estimated by measuring the concentration changes along HDO gradients in artificial ice in the range of densities ≈ 580–800 kg m⁻³. The measurements presented in that study are consistent with a close-off density for the diffusive fluxes equal to ≈ 800–805 kg m⁻³. However due to the lack of data covering the range of the lower densities we perform here a sensitivity test investigating the influence of the tortuosity uncertainty on the diffusion length value for type 1 and type 2 climate forcing regimes.

In the Iso-CFM it is possible to adjust the value of ρ_co. This in turn results in modifying the value of the scaling parameter b _τ in Eqn (19) such that 1 − b _τ(ρ_co/ρ_ice) = 0. We perform 2000 repetitions of a steady-state run using the HLD model where the ρ_co value is drawn from a normal distribution with a mean of 804.3 kg m⁻³ and a SD of 15 kg m⁻³. We ran the experiment for the type 1 and type 2 climate forcing regimes. The resulting diffusion length distributions (σ′₁₈) are presented in m of ice eq. (Fig. 8). Based on the results of the calculation the uncertainty (1-SD) on the diffusion length is ≈ ±0.0025 m for type 1 and ± 0.0034 m for the type 2 climate forcing (Fig. 8). Moreover, the influence of the varying close-off density and thus the tortuosity profile is increasingly more important with increasing depth. Based on the contour plot of Figure 3 the resulting uncertainty with respect to temperature is also similar for both climate forcing regimes and equal to ≈ 1 K. The test indicates the importance of the tortuosity profile and its impact on the diffusivity coefficient and diffusion length profiles.

Fig. 8. (a) Type 1 climate forcing. (b) Type 2 climate forcing. Close-off density/tortuosity sensitivity test. 2000 repetitions with ρ_co = 804.3 ± 15kg m⁻³. The HLD densification model is used for the test. The 1-standard deviation intervals of inverse tortuosity, the diffusivity and the diffusion length are represented by the black lines while the red lines represent the mean profiles of these quantities. The values of σ′₁₈ are given in m ice eq. and ‘b coefficient’ refers to the tortuosity scaling factor $b_\tau$ as given in Eqn (19) and Table 6.

Influence of the amplitude of the temperature seasonal cycle

To investigate the influence of the amplitude of the temperature seasonal cycle, we perform an experiment where we run all four implementations of the model for type 1 and type 2 conditions, for a range of values of the T _amp parameter in Eqn (32). We vary T _amp in the range [0, 14] K with a step of 1 K. For the case of type 1 conditions, we perform a spin-up run of 1000 years and a main run of 2000 years, while for the type 2 regime the spin-up and main runs are 500 years long each.

As seen in Figure 9, the influence of the T _amp parameter on the value of σ′_D, σ′₁₈ and σ′₁₇ is apparent, although minimal, while the value of the diffusion length increases with T _amp due to the nonlinear dependence of the diffusivity coefficient to temperature. The difference between diffusion lengths in the T _amp = 0 K and T _amp = 14 K runs is on the order of 5 × 10⁻⁴ m for all models and all diffusion lengths, translating to a temperature difference in the order of mK. Additionally, all four models behave similarly for both type 1 and type 2 conditions. For the GOU model we note that for the type 2 (higher temperature and accumulation) conditions, its deviation is higher than the other three models, a finding that is in line with the previous Sections ‘Steady-state climate forcing’ and ‘Diffusion length profiles’.

Fig. 9. Seasonal cycle amplitude test for type 1 (subplots a, b, c) and type 2 (subplots d, e, f) conditions. The Δσ′ notation represents the difference $\sigma '\left(T_{{\rm amp}} = T_{{\rm amp}}^k\right)- \sigma '\left(T_{{\rm amp}} = 0 \right)$ with k ∈ [0,  14 K] and T _amp the amplitude of the seasonal temperature signal.

It should be mentioned that the densification models tested here, even though typically used for ice core studies may not be the best choices for investigating the influence of the seasonal cycle of temperature. Preliminary results from densification studies using firn compaction data from South Pole and the EastGRIP sites, indicate that the four models tested in this study tend to underestimate the influence of the seasonal signal on the densification process. Further tests including other densification modelling approaches – some of them already present in Iso-CFM – can provide better answers.

Ramp experiments

We test the performance of the diffusion models for the case of a transient simulation where the climate forcing ramps from low to high temperature–accumulation conditions. The total time of the simulation is 10 000 years and it consists of 4000 years of steady-state conditions with T = 233.15 K and $A = 0.1 \,{\rm m\, a}^{-1}$ ice eq., a linear climatic transition to T = 248.15 K and $A = 0.2 \, {\rm m\, a}^{-1}$ ice eq. spanning 2000 years and steady-state conditions thereafter. The heat diffusion module of the model is enabled. The forcing of the experiment is illustrated in Figure 10a. In Figure 10b we present the close-off depth (ρ_co = 804.3 kg m⁻³) for all four models at every time step of the simulation and the results for σ′_D, σ′₁₈ and σ′₁₇ are given in Figure 11. In addition to the four model implementations, the analytical solution using the Herron and Langway steady-state implementation for σ′₁₇, σ′₁₈ and σ′_D is evaluated for the climate forcing of the experiment and illustrated as well.

Fig. 10. (a) Temperature and accumulation forcing for the transient simulation. (b) The close-off depth for the transient simulation and the age at the close-off depth modelled with the BAR model.

Fig. 11. Transient simulation results for the diffusion length at the close-off depth expressed in m of firn at the close-off density: (a) σ′_D, (b) σ′₁₈, (c) σ′₁₇. The black line represents the diffusion length value as calculated using the steady-state analytical expressions from Eqns (29) and (30).

As in the results of Section ‘Steady-state climate forcing’, the first 4000 years of the simulation $\lpar t_{\rm model} \in \lsqb -10\, 000\comma\; -6000\, {\rm years}\rsqb \rpar$ under steady-state climate, reveal the differences between the densification rate sensitivities of the different models to temperature and accumulation. An interesting aspect of the ramp experiment is the observation of the response time of the models. As the warming signal diffuses into the firn, one should expect a delay at the onset of the temperature and accumulation rise between the forcing and the σ′ signals. A difference of ≈ 400 years is observed between the onset of the forcing transition and the first upward change in the value of σ′₁₈. The magnitude of the delay does not vary between the different diffusion length signals (σ′_D, σ′₁₈ and σ′₁₇) with the delay itself being slightly shorter than the age of the close-off depth (≈500 years – Fig. 10b), pointing to the combined effect of the advection and diffusion acting more effectively in transmitting the warming signal in the firn column. As the climate forcing reaches steady-state conditions at t _model = −4000 years the response of the diffusion length signal is delayed by only ≈ 175 years, a number that reflects the shallower close-off depth for the higher temperature–accumulation conditions. A closer look into the firn temperature profile signal reveals a temperature gradient that persists for more than 2000 years after $t_{\rm model} = -4000 \, {\rm years}$. This diminishing temperature gradient results in a very slowly changing value of σ′₁₈ through this time interval, which does not settle until the end of the simulation time (t = 0). For a visualisation of this effect the reader is referred to the animation HLdynamic_ramp.mp4 included in the SOM.

Additionally, we present the effects of the temperature and accumulation rate change separately on both the diffusion length and the close-off depth. In Figure 12, the solid lines show the effect of temperature in increasing the diffusion length and shallowing the close-off depth as the accumulation is kept constant and the temperature increases. Similarly, the dash lines represent the experiments for which we isolate the impact of the accumulation rate change while the temperature is kept constant. It is apparent that the temperature change of 15 K has a greater impact on the value of the diffusion length compared to that of the doubling of the accumulation rate. We also observe that when the accumulation is doubled while the temperature is kept constant, the differences between the four models are far greater for the close-off values whereas the diffusion length results are comparable for the four model implementations.

Fig. 12. Individual effect of the temperature and the accumulation forcing on the diffusion length (a) and the close-off depth (b). Solid (dash) lines represent the experiments for which the accumulation (temperature) has been kept constant.

Warming pulse experiments

We investigate the models’ response to a warming pulse with a forcing that consists of an initial steady-state period, followed by a rapid warming and cooling phase. The pulse lasts 2000 years, has a temperature magnitude of 10 K (223.15 − 233.15 K) and an accumulation magnitude of $0.05 \, {\rm m\, a}^{-1}$ ice eq. $\lpar 0.05 - 0.1 \, {\rm m\, a}^{-1}\rpar$. In order to assess the evolution of σ′₁₈ under various warming and cooling rates, we introduce a characteristic time ψ that controls the rate of change for the temperature and accumulation forcing signals at the onset and end of the pulse. The forcing signal can be described as (Eqns 33 and 34):

(33)$$\eqalign{T_{{\rm forcing}} & = 233.15 \,\,\,\,\,-1000 \lt t_{{\rm model}} \le 0\cr T_{{\rm forcing}} & = 233.15 + 10 \left( 1 - \exp \left( - {t_{{\rm model}} - t_{on}\over \psi} \right) \right) \cr &\,\,\,\,\,0 \lt t_{{\rm model}} \le 2000\cr T_{{\rm forcing}} & = 233.15 + T_{{\rm off}} \exp \left( - {t_{{\rm model}} - t_{{\rm off}}\over \psi} \right) \cr &\,\,\,\,\,2000 \lt t_{{\rm model}} \le {10\, 000}}$$

and

(34)$$\eqalign{A_{{\rm forcing}}& = 0.05 \,\,\,\,\,-1000 \lt t_{{\rm model}} \le 0\cr A_{{\rm forcing}}& = 0.05 + 0.05 \left( 1 - \exp \left( - {t_{{\rm model}} - t_{on}\over \psi} \right) \right) \cr &\,\,\,\,\,0 \lt t_{{\rm model}} \le 2000\cr A_{{\rm forcing}} & = 0.05 + A_{{\rm off}} \exp \left( - {t_{{\rm model}} - t_{{\rm off}}\over \psi} \right)\cr & \,\,\,\,\,2000 \lt t_{{\rm model}} \le {10\, 000}}$$

We implement five forcing scenarios with ψ = 50, 75, 100, 200, 300, years. For all scenarios, $t_{{\rm on}} = 0 \, {\rm year}$, $t_{{\rm off}} = 2000 \, {\rm years}$, and T _off and A _off are determined by the value of T _forcing and A _forcing at $t_{{\rm model}} = 2000 \, {\rm years}$. The spin-up run is 4000 years long (−5000 < t _model < −1000, not shown in Figure 13) and the thickness of the firn column under consideration is 200 m. We use a time step of 1 year, ρ₀ = 350 kg m⁻³ and P = 0.7 Atm. We present the results of the warming pulse experiment in Figures 13a–d, and in the animation HLdynamic_RC_tau_100.mp4 in the SOM. Based on the outcome of the simulations, the following observations can be made.

Fig. 13. Results of pulse test for BAR, GOU, HLD and HLS models (subplots a–d), for $\psi = 50\comma\; \, 75\comma\; \, 100\comma\; \, 200\comma\; \, 300 \, {\rm years}$. The solid lines represent the temperature forcing while lines with bullet markers represent the diffusion length signal.

For all the simulations there is a time delay between the onsets of the warming signal and the diffusion length signal σ′₁₈. The delay is similar between all the models considered and approximately equal to the age of the close-off depth which for the warm phase climate conditions of the simulation is ≈400 years. We also observe that longer characteristic times ψ, result in an increasingly longer delay, reflecting the combined effect of the slower shallowing of the firn column as well as the slower penetration of heat in the firn due to the less steep temperature and accumulation gradients.

The test reveals the differences between the models’ dynamic behaviour under changing climate forcing at various rates. Overall, the models show a comparable behaviour across all the rates of change scenarios, although a notable difference can be observed in the initial fast decrease of the diffusion length value at the onset of the warming for the GOU and BAR models, in the order of 10⁻³ m. Similarly, the same models overshoot at the beginning of the cooling phase of the pulse, a behaviour that is followed within 100–200 years by the expected behaviour of a decreasing diffusion length during the cooling phase of the pulse.

This seemingly counter-intuitive result depends on the time constant ψ and is more profound for pulses of higher rapidity (low ψ values). Simulations where we have separated the temperature and accumulation effects on the diffusion length signal show that the BAR and particularly the GOU model result in a more rapid decrease in the σ′₁₈ signal when we only change the accumulation and keep the temperature constant. On the other hand, when we introduce pulses of temperature keeping the accumulation constant, the rate of change of the σ′₁₈ signal is very similar for all four models (see Figs 20 and 21 in the Appendix) This initial, accumulation-induced, rapid decrease in the σ′₁₈ signal overcomes the effect of temperature especially for the GOU model although its duration is short-lived and is followed by an increase in the σ′₁₈ signal, concomitant to the temperature increase.

A closer look at the warming phase of the experiment reveals that the σ′₁₈ value does not reach steady-state for any of the models and characteristic times considered. Later, as the temperature and the accumulation decreases at $t_{{\rm model}} = 2000\, {\rm years}$, the value of σ′₁₈ decreases with a long decay time. This decay is due to the residual heat in the firn column, which results in a temperature gradient that persists for $\approx \!6000\, {\rm years}$. This effect can be of importance for the study of the σ′₁₈ signal at millennial time scales and for the purpose of reconstructing past temperatures. The use of an analytical model approach based on the inversion of Eqns (29) and (30) for temperature and assuming steady-state will result in a less accurate temperature estimation over strong and abrupt climatic transitions.

Estimating signal attenuation from the δ¹⁸O record

In order to estimate the annual signal attenuation based on the high-resolution δ¹⁸O record, we make use of the analytical integration technique for spectral peaks described in Johnsen and Andersen (Reference Johnsen and Andersen1978). We calculate the power spectral density of 5 m sections of the δ¹⁸O signal using the maximum entropy method (MEM hereafter) (Burg, Reference Burg1975). For every 5 m section we interpolate the δ¹⁸O time series on the mean resolution of the section. Subsequently, using the autoregressive coefficients of the MEM prediction filter, we perform an analytical calculation of the position (in the frequency space) and the total area of the annual spectral peak. This in turn yields an estimate of the annual layer thickness λ_A (in other words the wavelength of the δ¹⁸O annual signal) and the mean magnitude of the annual component (hereafter Γ_A) in the δ¹⁸O signal for every 5 m depth interval. The value of Γ_A for the top depth interval is equal to $\Gamma _{{\rm A0}} = 2.5\hbox {\textperthousand } \permil$ and is used to obtain a data-based estimate of the signal attenuation as a function of depth.

We perform a steady-state 1000 years run in the Iso-CFM for every of the four model implementations considered here. Estimates of the density profiles are compared to the measured densities from the site (Fig. 15). The comparison reveals a profound mismatch between the GOU model and the measured density profile, with the model predicting faster densification rates throughout the full firn column. The HLD and HLS implementations show the best agreement whereas the BAR model deviates from the data for densities above 800 kg m⁻³. Using the modelled σ₁₈ signal we calculate the decay of the magnitude of the annual signal $\Gamma _{{\rm A}}/\Gamma _{{\rm A0}} = \exp \left(-2\pi ^2{\sigma ^2_{18}}/\lambda ^2_{\hbox{A}} \right)$ where λ_A is the spline-smoothed annual layer thickness based on the counted chronology.

Fig. 15. Measured density profile for the Site-A shallow core and CFM modelled profiles.

The results of the model and data-based spectral estimation of the annual signal magnitude decay are shown in Figure 16. We see that the HLD and HLS type models are the most accurate in predicting the attenuation of the annual signal. The observed and modelled annual signal magnitudes, decay exponentially with depth until the close-off depth after which they reach a constant value. The faster densification rates predicted by the GOU model yield a lower value for the diffusion length σ₁₈. As a result the GOU-inferred magnitude of the annual signal is higher than all the other models as well as the data-based magnitude. The slight discrepancy between the BAR model and the measured density at the bottom of the core results in a mismatch between the BAR-modelled and the data-based annual signal magnitude. This result shows how the densification of firn can have an effect on the diffusion length signal even past the close-off depth (≈ 70.9 m) and thus in parts of the firn where the diffusion process has ceased.

Fig. 16. Annual signal magnitude decay and annual layer thickness λ_A for the Site-A shallow core. We estimate independently the annual layer thickness λ_A from (1) the counted annual layers fitted with a smoothing spline and (2) by means of the spectral peak detection/integration technique in the MEM power spectrum (Johnsen and Andersen, Reference Johnsen and Andersen1978) at 5 m resolution sections. The black solid line in the bottom plot represents the magnitude of the annual signal as estimated using the spectral peak analytical integration method whereas the coloured lines are calculated using the Iso-CFM σ₁₈ profiles in combination with spline smoothed annual layer thickness λ_A profile based on the counted chronology.

Description of the ice core datasets

The ice core sites considered are Dome F, the Dome C and the EDML. For site characteristics and data references see Table 2. All the sections considered here contain data from depths below the close-off, though still relatively shallow compared to the full length of the cores. This has a threefold advantage: (a) the effect of solid ice diffusion is negligible, (b) the ice flow thinning is minimal and accurately constrained and (c) the accumulation rate for all three sites shows very small deviations compared to present for the time windows considered (Veres and others, Reference Veres2013).

Table 2. Ice core data sections and the corresponding drill site characteristics

A summary of the datasets can also be found in Holme and others (Reference Holme, Gkinis and Vinther2018).

Source of data: Svensson and others (Reference Svensson2015)^a, Gkinis (Reference Gkinis2011)^b, Oerter and others (Reference Oerter, Graf, Meyer and Wilhelms2004)^c.

The three datasets have each been measured with a different analytical technique. For the EDML and Dome C sections, two different Isotope Ratio Mass Spectrometry techniques were used (Oerter and others, Reference Oerter, Graf, Meyer and Wilhelms2004; Gkinis, Reference Gkinis2011). The Dome F section has been analysed with Cavity Ring Down Spectrometry using a continuous flow analysis approach (Gkinis and others, Reference Gkinis2011; Svensson and others, Reference Svensson2015). For all sections, δ¹⁸O and δD data are available, thus allowing for a comparison between the reconstructions of the two diffusion length signals. Moreover, two published temperature reconstruction studies (Stenni and others, Reference Stenni2010; Uemura and others, Reference Uemura2012) based on the isotope mixed cloud model (hereafter IMCM) (Ciais and Jouzel, Reference Ciais and Jouzel1994) are available, allowing for a comparison of our results with the more traditional water isotope thermometer method using the combined δD, D_excess signal to infer site–source temperature variations.

The two East Antarctic sites included in this test (Dome F and Dome C) are particularly interesting as they are characterised by a very similar annual accumulation signal while there is a difference of ~4 K in the present temperature signal. A strong accumulation intermittency is also observable with only a few events comprising the total annual accumulation. Although the prevalent conditions at the site are cold and dry, with typical occurrence of clear sky precipitation, a significant part of the snowfall takes place under rare episodes of warmer conditions (Fujita and Abe, Reference Fujita and Abe2006). These features can possibly pose a challenge to the classical isotope palaeothermometer (Fujita and Abe, Reference Fujita and Abe2006; Schlosser and others, Reference Schlosser2017) by introducing biases due to (a) the non-representative temperatures under which precipitation is formed (b) the occurrence of kinetic effects observed under such cold and dry conditions. On the other hand, the diffusion process that we consider here is not affected by such effects, as the diffusion length for a certain temperature is independent of the initial isotopic signal at the surface. Although the traditional isotopic thermometer ‘samples’ individual precipitation events, diffusion is continuously affected by the temperature and the accumulation signal at the surface and the firn column.

Methodology and results of the experiment

The diffusion length estimates used for the three sections considered are presented in Table 4 of Holme and others (Reference Holme, Gkinis and Vinther2018). These values are given in m of ice equivalent and have been corrected for the effects of ice diffusion and sampling mixing as well as for the ice flow thinning. For a detailed description on the diffusion length estimation from high-resolution isotope data and the necessary corrections involved, the reader is referred to the methods sections in Gkinis and others (Reference Gkinis, Simonsen, Buchardt, White and Vinther2014) and Holme and others (Reference Holme, Gkinis and Vinther2018).

We convert the ice equivalent diffusion length values from Table 4 in Holme and others (Reference Holme, Gkinis and Vinther2018) to m of firn equivalent at the depth of the close-off density ρ _co = 804.3 kg m⁻³ by scaling with the factor ρ _ice/ρ _co. This scaling yields the data-based estimate of the diffusion length $\hat {\sigma }'$. Subsequently, we estimate the root of the equation $\sigma '\lpar T\comma\; \, \rho = \rho _{{\rm co}}\rpar - \hat {\sigma }' = 0$ from which an absolute temperature estimate in K is obtained. For the root finding step we use Brent's algorithm (Brent, Reference Brent1973; Press and others, Reference Press, Teukolsky, Vetterling and Flannery2007) and require a tolerance for the estimated root equal to 0.01 K. Every calculation of σ′(T,  ρ = ρ _co) is based on a 2500 years run of the model. We use the uncertainties of the diffusion length estimation from Table 4 in Holme and others (Reference Holme, Gkinis and Vinther2018) and generate normal distributions for the value of $\hat {\sigma }'$, which in turn yield a distribution of temperatures for each ice core site. The size of the distributions is equal to 500. For all our calculations the fractionation factors from Majoube (Reference Majoube1971) and Merlivat and Nief (Reference Merlivat and Nief1967) are used.

The mean and SDs of $\hat {\sigma }'$ and the estimated temperatures for all sites and all models are shown in Tables 3 and 4 and their relative distributions are presented in Figure 17. Overall, the results of the experiment indicate that the combined uncertainty of the four models is in the order of 0.5 K for all sites considered. The two isotopes yield very similar results (within 1-SD) for the Dome C and EDML sites, whereas the Dome F reconstruction indicates a slightly larger discrepancy between the $\hat {\sigma }'_{{\rm D}}$ and the $\hat {\sigma }'_{{\rm 18}}$ based temperature estimates (marginally over 2-SDs).

Fig. 17. Results of the temperature estimation test for three Antarctic sites. For every core and every isotopologue we present the starting distribution for the diffusion length signal as found in Holme and others (Reference Holme, Gkinis and Vinther2018) and scaled to its close-off density value σ′₁₈ (m of firn eq.) as well as the resulting temperature distributions for all four models. Dot-dash lines represent the temperature estimate based on the IMCM and dot lines represent the mean temperature estimate based on both the $\hat {\sigma }'_{{\rm D}}$ and the $\hat {\sigma }'_{{\rm 18}}$ reconstructions.

Table 3. Results of the ice core data experiments for δD – $\hat {\sigma }_{{\rm D}}$

Table 4. Results of the ice core data experiments for δ¹⁸O – $\hat {\sigma }_{{\rm 18}}$

In Table 5 we present the comparison of our results with the present temperature as well as the temperature estimates for the time windows under consideration based on the IMCM. The youngest age of the Dome F and Dome C temperature reconstructions is 100 years B.P. (Stenni and others, Reference Stenni2010; Uemura and others, Reference Uemura2012) and is used as equivalent to present conditions for our comparison. The EDML reconstruction on the other hand has a youngest age of 1300 years B.P. (Stenni and others, Reference Stenni2010), so in order to compare to present, we use the regional reconstruction from PAGES 2k Consortium (2013) and account for a mean temperature anomaly of 0.43 K between present and the reference period 1.2–2 ka B.P. in Stenni and others (Reference Stenni2010). The comparison is presented in Figure 17 where the dot and dot-dash vertical lines represent the diffusion and the IMCM temperature estimates respectively. The two techniques agree within 1-SD (we refer to the SD of the diffusion reconstruction. The IMCM estimates do not have an adjoint uncertainty interval) for the Dome F and Dome C sites. On the other hand, the EDML reconstructions differ by 1.1 K, a result that lies outside of the 3-SD interval of the diffusion estimate and which we comment further on, in Section ‘Antarctic sites temperature reconstruction – comparison with the isotope mixed cloud model’.

Table 5. Comparison between the diffusion and the IMCM temperature reconstructions

For Dome F the IMCM temperature reconstruction is from Uemura and others (Reference Uemura2012)^a, for Dome C we have used Stenni and others (Reference Stenni2010)^b and for EDML we combined the reconstruction from Stenni and others (Reference Stenni2010)^b and the PAGES 2k regional reconstruction (PAGES 2k Consortium, 2013). The column T _{diffusion–all} contains the combined T _D–all, T _18–all diffusion-based temperature estimates.