2.1.1 Step 1: prior input

The starting point of the optimization procedure is the first guess. To construct the first-guess temperature input T_g,0(t), a constant temperature of −29.6 ^∘C is used for the complete Holocene section, which corresponds to the last value of the temperature spin-up (Fig. 2b).

2.1.2 Step 2: Monte Carlo type input generator – generating long-term solutions

During the second step of the optimization, the prior temperature input T_g,0(t) from step 1 is iteratively (j) changed following a Monte Carlo approach. The basic idea of the Monte Carlo approach is to generate smooth temperature inputs T_mc,j(t) by superimposing low-pass-filtered values P_j of uniformly distributed random values P_r,j on the prior input $T_{\mathrm{mc},j-\mathrm{1}}$. Then, the new input is fed to the firn model and the mismatch $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},j}}$ (with X≡δ¹⁵N_mc,j) between the modelled δ¹⁵N_mc,j (here X_mod), calculated from the model output, and the synthetic δ¹⁵N_syn (here X_target) is computed for every time step (i) of the target data δ¹⁵N_syn according to

(Note: if not otherwise stated, all mismatches in this study labelled with “D” are calculated similar to Eq. 1.)

Table 2Summary for the Monte Carlo approach: mismatch D_g,0 between the modelled δ¹⁵N (or temperature) values using the first-guess input and the synthetic δ¹⁵N (or temperature) target for each scenario. D_mc is the mismatch between the modelled δ¹⁵N (or temperature) using the final Monte Carlo temperature solution and the synthetic δ¹⁵N (or temperature) target for each scenario.

Download Print Version | Download XLSX

$D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc}}}$ serves as the criterion which is minimized during the optimization in step 2. If the mismatch $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},j}}$ decreases compared to the prior input ($T_{\mathrm{mc},j-\mathrm{1}}$, $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},j-\mathrm{1}}}$), the new input is saved and used as the new guess ($T_{\mathrm{g},j}=T_{\mathrm{mc},j}$). This procedure is repeated until convergence is achieved, leading to the final long-term temperature T_mc,fin(t). Table 2 lists the number of improvements and iterations performed for the different synthetic datasets.

The perturbation of the current guess T_g,j is conducted in the following way: let ${\mathit{T}}_{\mathrm{g},\mathrm{0}}=T_{\mathrm{g},\mathrm{0}}\left(t\right)$ be the vector containing the prior temperature input. A second vector (P_r,1) with the same number of elements n_mc as T_g,0 is generated, containing n_mc uniformly distributed random numbers within the limits of an also randomly (equally distributed) chosen standard deviation (SD). SD is chosen from a range of 0.05–0.50, which means that the maximum allowed perturbation of a single temperature value T(t₀) is in a range of ±5 to ±50 %. Creating the synthetic frequencies, P_r,1 is low-pass filtered using cubic-spline filtering (Enting, 1987) with an equally distributed random cut-off period (COP) in the range of 500 to 2000 years generating the vector P₁. Hereby the low-pass filtering of P_r,1 reduces the amplitudes of the perturbation of T_g,0. The new surface temperature input T_mc,1 is calculated from P₁ according to

The superscript “T” stands for transposed and $\widehat{\mathrm{1}}$ is the n by 1 matrix of ones.

This approach provides a high potential for parallel computing. In this study, an eight-core computer was used, generating and running eight different inputs of T_mc simultaneously, minimizing the time to find an improved solution. For example, during the 706 iterations for scenario S2, about 5600 different inputs were created and tested, leading to 351 improvements (see Table 2). Since it is possible to find more than one improvement per iteration step due to the parallelization on eight CPUs, the solution giving the minimal misfit $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},j}}$ is chosen as new first-guess for the next iteration step. This leads to a decrease of the used improvements for the optimization (e.g. for S2, 172 of the 351 improvements were used). Additionally, a first gas-age scale (Δage_mc,fin(t)) is extracted from the model using the last improved conditions, which will then be used in step 3.

2.1.3 Step 3: adding short-term (high-frequency) information

In step 3, the missing short-term temperature history providing a suitable fit between modelled and synthetic δ¹⁵N data is directly extracted from the pointwise mismatch $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$, between the modelled δ¹⁵N_mc,fin(t) obtained in step 2 and the synthetic δ¹⁵N_syn target. Note that for a real reconstruction, this mismatch is calculated using the measured δ¹⁵N_meas dataset instead of the synthetic one. $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$ can be interpreted in first order as the detrended high-frequency signal of the synthetic δ¹⁵N_syn target. $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$ is transferred to the gas-age scale using Δage_mc,fin(t) provided by the firn-model output for the smooth temperature input T_mc,fin(t). This is needed to insure synchroneity between the high-frequency temperature variations ΔT(t) extracted from the mismatch $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$ on the ice-age scale and the smooth temperature solution T_mc,fin(t). Additionally, the signal is shifted by about 10 years towards modern values to account for gas diffusion from the surface to the lock-in depth (Schwander et al., 1993), which is not yet implemented in the firn model. This is necessary for adding the calculated short-term temperature changes ΔT(t) to the smooth signal T_mc,fin(t). The ΔT values are calculated according to Eq. (3):

using the thermal-diffusion sensitivity ${\mathrm{\Omega }}_{{\mathrm{N}}_{\mathrm{2}},i}$ for nitrogen-isotope fractionation from Grachev and Severinghaus (2003):

${\stackrel{\mathrm{‾}}{T}}_i$ is the mean firn temperature in Kelvin which is calculated by the firn model for each time point i. To reconstruct the final (high-frequency) temperature input (T_hf(t)), the extracted short-term temperature signal ΔT(t) is simply added to the long-term temperature input T_mc,fin(t):

2.1.4 Step 4: final correction of the surface temperature solution

For a further improvement of the remaining δ¹⁵N and resulting surface-temperature misfits ($D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$, $D_{T_{\mathrm{hf}}}\left(t\right)$), it is important to find a correction method that contains information that is also available when using measured data. The benefit of the synthetic data study is that several later-unknown quantities can be calculated and used for improving the reconstruction approach (see Sects. 3 and 4). For instance, it is possible to split the synthetic δ¹⁵N_syn data in the gravitational and thermodiffusion parts or to use the temperature misfit, which is unknown in reality. The idea underlying the correction algorithm explained hereafter is that the remaining misfits of δ¹⁵N ($D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$) and temperature ($D_{T_{\mathrm{hf}}}\left(t\right)$) are connected to the Monte Carlo (step 2) and high-frequency part (step 3) of the reconstruction algorithm. In the present inversion framework, it is not possible to find a long-term solution δ¹⁵N_mc,fin (or T_mc,fin) which exactly passes through the δ¹⁵N_syn (or T_syn) target in the middle of the variance in all parts of the time series. This leads to a slight over- or underestimation of δ¹⁵N_mc,fin(t) and their corresponding temperature values T_mc,fin(t). For example, a slightly too low (or too high) smooth temperature estimate T_mc,fin leads to a small increase (or decrease) of the firn-column height, creating a wrong gravitational background signal in δ¹⁵N_mc,fin on a later point in time (because the firn column needs some time to react). An additional error in the thermal-diffusion signal is also created due to the high-frequency part of the reconstruction (step 3), because the high-frequency information is directly extracted from the deviation of the synthetic target δ¹⁵N_syn(t) and the modelled δ¹⁵N_mc,fin(t) from the final long-term solution T_mc,fin(t) of the Monte Carlo part. Therefore, this error is transferred into the next step of the reconstruction and partly creates the remaining deviations.

Figure 3Scenario S1: (a) cross-correlation function (xcf) between IF(t) and the remaining mismatch in δ¹⁵N ($D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$), after the high-frequency part shows two extrema: the maximum correlation (max xcf) and the minimum correlation (min xcf). (b) Cross-correlation function (xcf) between IF(t) and the remaining mismatch in temperature ($D_{T_{\mathrm{hf}}}\left(t\right)$), after the high-frequency part, shows two extrema: the maximum correlation (max xcf) and the minimum correlation (min xcf). (c) Inverting of panel (a) in x (lag) and y (correlation coefficient) direction. (d) Comparison between panels (a) and (b). (e) Comparison between panels (a) and (c). The temperature-correction values are calculated from the linear dependency between IF(t) and $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$. After shifting IF(t) to max xcf (lag max) and to min xcf (lag min), Δδ¹⁵N_max(t) and Δδ¹⁵N_min(t) are calculated. Next, Δδ¹⁵N_max(t) and Δδ¹⁵N_min(t) are inverted. That means, for Δδ¹⁵N_max(t), the values are shifted back (−lag max) and shifted further to lag min. After inverting, Δδ¹⁵N_max(t) and Δδ¹⁵N_min(t) are summed up component-wise to calculate Δδ¹⁵N_cv(t). Using Δδ¹⁵N_cv(t) in Eq. (3) leads to the temperature-correction values which are added to the temperature T_hf.

Download

To investigate this problem, the deviations $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$ of the synthetic target data δ¹⁵N_syn to δ¹⁵N_mc,fin of the Monte Carlo part are numerically integrated over a time window of 200 years (Sect. 4, Supplement Sect. S3), and thereafter the window is shifted from past to future in 1-year steps resulting in a time series called IF(t). IF(t) equals a 200-year running mean of $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin}}}\left(t\right)$. For t, the middle position of the window is allocated. The time evolution of IF(t) is a measure for the deviation of the long-term solution δ¹⁵N_mc,fin(t) (or T_mc,fin(t)) from the perfect middle passage through the target data δ¹⁵N_syn(t) (or T_syn(t)) and for the slight over- and underestimation of the resulting temperature.

Next, the sample cross-correlation function (xcf) (Box et al., 1994) is applied to IF(t) and the remaining misfits $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$ of δ¹⁵N after the high-frequency part. The xcf shows two extrema (Fig. 3a), a maximum (xcf_max) and a minimum (xcf_min) at two certain lags (lag${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ at xcf${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ and lag${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ at xcf${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$). Now, the same analysis is conducted for IF(t) versus the temperature mismatch $D_{T_{\mathrm{hf}}}\left(t\right)$ (Fig. 3b), which shows an equal behaviour (two extrema, lag${}_{max,T}$ at xcf${}_{max,T}$ and lag${}_{min,T}$ at xcf${}_{min,T}$). Comparing the two cross correlations shows that lag${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ equals the negative lag${}_{min,T}$ and lag${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ corresponds to the negative lag${}_{max,T}$ (Fig. 3d, e). The idea for the correction is that the extrema in the cross-correlation IF(t) versus $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$ with the positive lag (positive means here that $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$ has to be shifted to past values relative to IF(t)) creates the misfit of temperature $D_{T_{\mathrm{hf}}}\left(t\right)$ on the negative lag (modern direction) of IF(t) versus $D_{T_{\mathrm{hf}}}\left(t\right)$, and vice versa. So, IF(t) yields information about the cause and allows us to correct this effect between the remaining mismatches, $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$ and $D_{T_{\mathrm{hf}}}\left(t\right)$, over the whole time series. The lags are not sharp signals, due to the fact that (i) the cross correlations are conducted over the whole analysed record, leading to an averaging of this cause-and-effect relationship and that (ii) IF(t) is a smoothed quantity itself. The correction of the reconstructed temperature after the high-frequency part is conducted in the following way: from the two linear relationships between IF(t) and $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$ at the two lags (lag${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ at xcf${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$, lag${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ at xcf${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$), two sets of δ¹⁵N correction values (Δδ¹⁵N_max(t) from xcf${}_{max,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ and Δδ¹⁵N_min(t) from xcf${}_{min,{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$) are calculated. Then, the lags are inverted (Fig. 3c, e), shifting the two sets of the δ¹⁵N correction values to the attributed lags of the cross correlation between IF(t) and $D_{T_{\mathrm{hf}}}\left(t\right)$ (e.g. Δδ¹⁵N_min(t) to lag from xcf${}_{max,T}$ from the cross correlation between IF(t) and $D_{T_{\mathrm{hf}}}\left(t\right)$), therefore changing the time assignments of Δδ¹⁵N_min(t) and Δδ¹⁵N_max(t) to Δδ¹⁵N${}_{min}(t+{\mathrm{lag}}_{max,T}$) and Δδ¹⁵N${}_{max}(t+{\mathrm{lag}}_{min,T}$). Now, the Δδ¹⁵N_max(t) and Δδ¹⁵N_min(t) are summed up component-wise, leading to the time series Δδ¹⁵N_cv(t). From Eq. (3) with Δδ¹⁵N_cv,i instead of $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc},\mathrm{fin},i}}$, the corresponding temperature correction values are calculated and added to the high-frequency temperature solution T_hf(t), giving the corrected temperature T_corr(t). Finally, T_corr(t) is used to run the firn model to calculate the corrected δ¹⁵N_corr(t) time series. This cause-and-effect relationship found in the cross correlations between IF(t) and $D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{hf}}}\left(t\right)$, and IF(t) and $D_{T_{\mathrm{hf}}}\left(t\right)$, is exemplarily shown in Fig. 3 for scenario S1 and was found for all eight synthetic scenarios. The derived correction algorithm leads to a further reduction of the mismatches of about 40 % in δ¹⁵N and temperature (see Sect. 3.2).

2.3.1 Timescale

For the entire study, the GICC05 chronology is used (Rasmussen et al., 2014; Seierstad et al., 2014). During the whole reconstruction procedure, the two input time series (surface temperature and accumulation rate) are split into two parts. The first part ranges from 20 to 10 520 years b2k (called the “Holocene section”) and the second one from 10 520 to 35 000 years b2k (“spin-up section”). The entire accumulation-rate input, as well as the spin-up section of the surface-temperature input, remains unchanged during the reconstruction procedure.

2.3.2 Accumulation-rate data

Besides surface temperatures, accumulation-rate data are needed to drive the firn model. In this study, we use the original accumulation rates, reconstructed in Cuffey and Clow (1997), produced using an ice-flow model adapted to the Greenland Ice Sheet Project Two (GISP2) location but adapted to the GICC05 chronology (Rasmussen et al., 2008; Seierstad et al., 2014). A detailed description of the adaption procedure can be found in Sect. S1 of the Supplement. The raw accumulation-rate data for the main part of the spin-up section (12 000 to 35 000 years b2k) are linearly interpolated to a 20-year grid and low-pass filtered with a 200-year COP using cubic-spline filtering (Enting, 1987). For the Holocene section (20–10 520 years b2k) and the transition part between Holocene and spin-up section (10 520 to 12 000 years b2k), the raw accumulation-rate data are linearly interpolated to a 1-year grid to obtain equidistant integer point-to-point distances which are necessary for the reconstruction, and to preserve as much information as possible for this time period (Fig. 2a). Except for these technical adjustments, the accumulation-rate input remains unmodified, assuming high reliability of these data during the Holocene. The accumulation data were reconstructed using annual layer counting, and a thinning model which should lead to maximum relative uncertainty of 10 % for the first 1500 of the 3000 m ice core (Cuffey and Clow, 1997). From the three accumulation-rate scenarios reconstructed in Cuffey and Clow (1997) and adapted here to the GICC05 chronology, the intermediate one is chosen (red curves in Fig. S1 in the Supplement). Since the differences between the scenarios are not important for the evaluation of the reconstruction approach, they are not taken into account for this study.

Additionally, two sensitivity experiments were conducted (see Sect. S2 in the Supplement) in order to investigate (i) the influence of low-pass filtering of the high-resolution accumulation rates on the model outputs and (ii) the possible contribution of the accumulation-rate variability to the δ¹⁵N data during the Holocene. The first experiment shows that filtering the accumulation rates with cut-off periods in the range of 20 to 500 years has nearly no influence on the modelled δ¹⁵N or lock-in depth as long as the major trends are being conserved. The second experiment leads to the finding that the accumulation-rate variability explains about 12 to 30 % of δ¹⁵N variability. A total of 30 % corresponds to the 8.2 kyr event and 12 % to the mean of the whole Holocene period including the 8.2 kyr event. Hence, the influence of accumulation changes, excluding the extreme 8.2 kyr event, is generally below 10 % during most parts of the Holocene.

2.3.3 δ¹⁸O_ice data

Oxygen-isotope data from the GISP2 ice-core-water samples measured at the University of Washington's Quaternary Isotope Laboratory are used to construct the surface-temperature input of the model spin-up (12 to 35 kyr b2k, Grootes et al., 1993; Grootes and Stuiver, 1997; Meese et al., 1994; Steig et al., 1994; Stuiver et al., 1995; data availability: Grootes and Stuiver, 1999). The raw δ¹⁸O_ice data are filtered and interpolated in the same way as the accumulation-rate data for the spin-up part.

2.3.4 Surface-temperature spin-up

The surface-temperature history of the spin-up section (Fig. 2b) is obtained by calibrating the filtered and interpolated δ¹⁸O_ice data (Eq. 13) using the values for the temperature sensitivity α_18O and offset β found by Kindler et al. (2014) for the North Greenland Ice Core Project (NGRIP) ice core assuming a linear relationship of δ¹⁸O_ice with temperature.

The values 35.2 ‰ and −31.4 ^∘C are modern-time parameters for the GISP2 site (Grootes and Stuiver, 1997; Schwander et al., 1997). The spin-up is needed to bring the firn model to a well-defined starting condition that takes possible memory effects (influence of earlier conditions) of firn states into account.

Figure 4(I) S1–S5: synthetic smooth temperature scenarios for the construction of the target-temperature data. (II) S1–S5: synthetic target surface-temperature scenarios. (III) S1–S5: corresponding synthetic δ¹⁵N target time series.

Download

Figure 5(I) Synthetic target surface-temperature scenarios H1–H3. (II) Corresponding synthetic δ¹⁵N target time series H1–H3. (III) GISP2 δ¹⁵N measured data (Kobashi et al., 2008b).

Download

2.3.5 Generating synthetic target data

In order to develop and evaluate the presented algorithm, eight temperature scenarios were constructed and used to model synthetic δ¹⁵N data, which serve later as targets for the reconstruction. From these eight synthetic surface-temperature and related δ¹⁵N scenarios (S1–S5 and H1–H3), three datasets (later called Holocene-like scenarios H1–H3) were constructed in such a way that the resulting δ¹⁵N time series are very close to the δ¹⁵N values measured by Kobashi et al. (2008b) in terms of variability (amplitudes) and frequency (data resolution) of the GISP2 nitrogen-isotope data (Figs. 4, 5).

The synthetic surface-temperature scenarios S1–S5 are created by generating a long-term temperature time series (T_syn,smooth) analogous to the Monte Carlo part of the reconstruction procedure for only one iteration step (see Sect. 2.1). The values for the cut-off period used for the filtering of the random values, and the SD values (standard deviation of the random values; see Sect. 2.1) for the first five scenarios can be found in Table 3. The long-term temperatures (Fig. 4I) are calculated on a 20-year grid, which is nearly similar to the time resolution of the GISP2 δ¹⁵N measurement values of about 17 years (Kobashi et al., 2008b). For the Holocene-like scenarios, the smooth temperature time series were generated from the temperature reconstruction for the GISP2 δ¹⁵N data (not shown here). The final Holocene surface-temperature solution was filtered with a 100-year cut-off to obtain the long-term temperature scenario.

Following this, high-frequency information is added to the long-term temperature histories. A set of normally distributed random numbers with a zero mean and a standard deviation (1σ) of 1 K for scenarios S1–S5 and 0.3 K for Holocene-like scenarios H1–H3 is generated on the same 20-year grid and added up to the long-term temperature time series. Finally, the resulting synthetic target-temperature scenarios (Figs. 4II, 5I) are linearly interpolated to a 1-year grid.

Table 3Cut-off periods (COPs) and SD values used for creating the smooth synthetic temperature scenarios according to the Monte Carlo approach.

Download Print Version | Download XLSX

Figure 6Evolution of the mean misfit ($D_{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{mc}}}$) of the modelled δ¹⁵N_mc versus synthetic target δ¹⁵N_syn as function of the number of iterations (j) for the Monte Carlo approach for all synthetic target scenarios.

Download

These synthetic temperatures are combined with the spin-up temperature and are used together with the accumulation-rate input to feed the firn model. From the model output, the synthetic δ¹⁵N targets are calculated according to Sect. 2.2. The firn-model output provides ice-age as well as gas-age information. The final synthetic δ¹⁵N target time series (δ¹⁵N_syn) are set intentionally on the ice-age scale to mirror measured data, because no prior information is available for the gas-age – ice-age difference (Δage) for ice-core data.

Figure 7(a–c) First-guess (g,0) versus Monte Carlo (mc,fin) δ¹⁵N solution for the scenario S5: (a) synthetic target δ¹⁵N_syn (black dotted line), modelled δ¹⁵N time series for the first-guess input (blue line) and Monte Carlo solution (red line). (b) Histogram shows the pointwise mismatches of $D_{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ for the first-guess solution (blue) and the Monte Carlo solution (red) versus the synthetic target. (c) Time series of the pointwise mismatches of $D_{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}$ for the first-guess solution (blue) and the Monte Carlo solution (red) versus the synthetic target. (d–f) First-guess versus Monte Carlo surface-temperature solution T_surf for the scenario S5: (d) synthetic surface-temperature target T_syn (black dotted line), first-guess temperature input (blue line) and Monte Carlo solution (red line). (e) Histogram shows the pointwise temperature mismatches D_Tfor the first-guess solution (blue) and the Monte Carlo solution (red) versus the synthetic surface-temperature target. (f) Time series of the pointwise temperature mismatches D_T for the first-guess solution (blue) and the Monte Carlo solution (red) versus the synthetic surface-temperature target.

Download

4.3.1 Benefits of the novel gas-isotope fitting approach

In addition to the fitting of δ¹⁵N data, the algorithm is able to fit δ⁴⁰Ar and δ¹⁵N_excess data as well using the same basic concepts (Fig. 10). Here, the δ⁴⁰Ar and δ¹⁵N_excess data from Kobashi et al. (2008b) were used as the fitting targets. We reach final mismatches (2σ) of 4.0 permeg for δ⁴⁰Ar/4 and 3.7 permeg for δ¹⁵N_excess, which are for both quantities below the analytical measurement uncertainty of 4.0 to 9.0 permeg for δ⁴⁰Ar/4 and 5.0 to 9.8 permeg for δ¹⁵N_excess measured data (Kobashi et al., 2008b).

The automated inversion of different gas-isotope quantities (δ¹⁵N, δ⁴⁰Ar, δ¹⁵N_excess) provides a unique opportunity to study the differences in the gained solutions using different targets and to improve our knowledge about the uncertainties of gas-isotope-based temperature reconstructions using a single firn model. Next, the presented algorithm is not dependent on the firn model, which leads to the implication that the algorithm can be coupled to different firn models describing firn physics in different ways. Furthermore, an automated reconstruction algorithm avoiding manual manipulation and leading to reproducible solutions makes it possible for the first time, to study and learn from the differences between solutions matching different targets. Finally, differences obtained by applying different firn physics (densification equations, convective zone, etc.) but the very same inversion algorithm may help to assess firn-model shortcomings, resulting in more robust uncertainty estimates than it was ever possible before.

In this publication, we show the functionality and the basic concepts of the automated inversion algorithm using well-known synthetic δ¹⁵N fitting targets. In this “perfect-world scenario”, the forward problem, converting surface temperature to δ¹⁵N, as well as the inverse problem, converting δ¹⁵N to surface temperature, are completely described by the used firn model. Consequently, all sources of signal noise are ignored. For the later use of the algorithm on δ¹⁵N, δ⁴⁰Ar or δ¹⁵N_excess measured data, this will not be the case anymore due to different sources of signal noise in the used measured data. As a result, differences between temperature solutions obtained from individual targets (δ¹⁵N, δ⁴⁰Ar, δ¹⁵N_excess) will become obvious. These differences will allow us to quantify the uncertainties associated with different unconstrained processes. Next, we will list and discuss potential sources of uncertainties and try to provide suggestions for their handling and quantification in our approach.

4.3.2 Measurement uncertainty and firn heterogeneity (centimetre-scale variability)

Many studies have investigated the influence of firn heterogeneity (or density fluctuations) on measurements of air compounds and quantities (e.g. δ¹⁵N, δ⁴⁰Ar, CH₄, CO₂, O₂ ∕ N₂ ratio, air content) extracted from ice cores resulting in centimetre-scale variability and leading to additional noise on the measured data (e.g. Capron et al., 2010; Etheridge et al., 1992; Fourteau et al., 2017; Fujita et al., 2009; Hörhold et al., 2011; Huber and Leuenberger, 2004; Rhodes et al., 2013, 2016). Using discrete measurement techniques instead of continuous sampling methods makes it difficult to quantify these effects. However, during discrete analyses of ice-core air data, it is common to measure replicates for given depths, from which the measurement uncertainties of the gas-isotope data are calculated using pooled standard deviation (Hedges and Olkin, 1985). Often, it is not possible to take real replicates (same depth), and instead the replicates are taken from nearby depths. Hence, any potential centimetre-scale variability is to some degree already included in the measurement uncertainty, because each measurement point represents the average over a few centimetres of ice. This is especially the case for low-accumulation sites or glacial ice samples for which the vertical length of a sample (e.g. 10–25 cm long for the glacial part of the NGRIP ice core; Kindler et al., 2014) covers the equivalent of 20 to 50 years of ice at approximately 35 kyr b2k. Increasing the depth resolution of the samples would increase our knowledge of centimetre-scale variability, for, e.g. identifying anomalous entrapped gas layers that could have been rapidly isolated from the surface due to an overlying high-density layer (e.g. Rosen et al., 2014). As this variability is likely due to heterogeneity in the density profile, modelling such heterogeneities (if possible at all) may not help to better reconstruct a meaningful temperature history but rather to reproduce the source of noise. This means that the potential centimetre-scale variability, in many cases, is already incorporated in the analytical noise obtained from gas-isotope measurements, due to analytical techniques themselves. Assuming the measurement uncertainty as Gaussian distributed, it is easy to incorporate this source of uncertainty in the inverse-modelling approach presented here. This will increase the uncertainty of the temperature according to Eq. (3).The same equation can also be used for the calculation of the uncertainty in temperature related to measurement uncertainty in general.

To answer the pertinent question of how to better extract a meaningful temperature history from a noisy ice-core record, an excellent – but costly – solution is of course to use multiple ice cores. For example, a δ¹⁵N-based temperature reconstruction from the combination of data from the GISP2 ice core with the “sister ice core” GRIP drilled 30 km apart is likely one of the best ways to overcome potential centimetre-scale variability. A comparison of ice cores that were drilled even closer might be even more advantageous.

4.3.3 Smoothing effects due to gas diffusion and trapping

It is known that gas-diffusion and trapping processes in the firn can smooth out fast signals and result in a damping of the amplitudes of gas-isotope signals (e.g. Grachev and Severinghaus, 2005; Spahni, 2003). The duration of gas diffusion from the top of the diffusive column to the bottom where the air is closed off in bubbles is for Holocene conditions in Greenland approximately of the order of 10 years (Schwander et al., 1997), whereas the data resolution of the synthetic targets was set to 20 years to mimic the measurement data from Kobashi et al. (2008b) with a mean data resolution of about 17 years (see Sect. 2.3: “Generating synthetic target data”). In the study of Kindler et al. (2014), it was shown that a glacial Greenland lock-in depth leads to a damping of the δ¹⁵N signal of about 30 % for a 10 K temperature rise in 20 years. We further assume that the smoothing according to the lock-in process is negligible for Greenland Holocene conditions according to the much smaller amplitude signals and shallower LID. Yet, for glacial conditions, it requires attention.

4.3.4 Accumulation-rate uncertainties

For the synthetic data study presented in this paper, it is assumed that the used accumulation-rate data are well known with zero uncertainty. This simplification is used to show the functionality and basic concepts of the presented fitting algorithm in every detail on well-known δ¹⁵N and temperature targets and to focus on the final uncertainties originating from the presented fitting algorithm itself. For the later reconstruction using measured gas-isotope data together with the published accumulation-rate scenarios shown in Sect. S1 of the Supplement, this will not be the case anymore. Uncertainties in layer counting and corrections for ice thinning will lead to a fundamental uncertainty. Especially in the early Holocene, this can easily exceed 10 %. As the accumulation-rate data are used to run the firn model, all potential accumulation uncertainties are in part incorporated into the temperature reconstruction. On the other hand, as we discussed in Sect. S2 of the Supplement, the accumulation-rate variability has a minor impact compared to the input temperature on the variability of δ¹⁵N data in the Holocene. The influence of these quantities, accumulation rate or temperature, on the temperature reconstruction is not equal; during the Holocene, accumulation-rate variability explains about 12 to 30 % of δ¹⁵N variability. A total of 30 % corresponds to the 8.2 kyr event and 12 % to the mean of the whole Holocene period including the 8.2 kyr event. Hence, the influence of accumulation changes, excluding the extreme 8.2 kyr event, is generally below 10 % during the Holocene. If the accumulation is assumed to be completely correct, then the missing part will be assigned to temperature variations. Nevertheless, for the fitting of the Holocene measurement data, we will use all three accumulation-rate scenarios as shown in Sect. S1 of the Supplement. The difference in the reconstructed temperatures arising from the differences of these three scenarios will be used for the uncertainty calculation as well and is most likely higher than the uncertainty arising from uncertainties due to the process of producing the accumulation-rate data and from the conversion of the accumulation-rate data to the GICC05 timescale.

4.3.5 Convective zone variability

Many studies have shown the existence of a well-mixed zone at the top of the diffusive firn column, called the convective zone (CZ). The CZ is formed by strong katabatic winds and pressure gradients between the surface and the firn (e.g. Kawamura et al., 2006, 2013; Severinghaus et al., 2010). The existence of a CZ changes the gravitational background signal in δ¹⁵N and δ⁴⁰Ar as it reduces the diffusive-column height. The presented fitting algorithm was used together with the two most frequently used firn models for temperature reconstructions based on stable isotopes of air, the Schwander et al. (1997) model which has no CZ built in (or better – a constant CZ of 0 m) and the Goujon firn model (Goujon et al., 2003) (which assumes a constant convective zone over time, that can easily be set in the code). This difference between the two firn models only changes significantly the absolute temperature rather than the temperature anomalies as it was shown by other studies (e.g. Guillevic et al., 2013, Fig. 3). In the presented work, we show the results using the model from Schwander et al. (1997), because the differences between the obtained solutions using the two models are negligible, except a constant temperature offset. Also noteworthy is that there is no firn model at the moment which uses a dynamically changing CZ. Indeed, this should be investigated but requires additional intense work. Additionally, the knowledge of the time evolution of CZ changes for time periods of millennia to several hundreds of millennia (in frequency and magnitude) is too poor to estimate the influence of this quantity on the reconstruction. In principle, it is possible to cancel out the influence of a potentially changing CZ by using δ¹⁵N_excess data for temperature reconstruction, as due to the subtraction of δ⁴⁰Ar/4 from δ¹⁵N the gravitational term of the signals is eliminated. From that point of view, it will be interesting to compare temperature solutions gained from δ¹⁵N_excess fitting with the solutions based on δ¹⁵N or δ⁴⁰Ar alone. This can offer a useful tool for quantifying the magnitude and frequency of CZ changes in the time interval of interest.

It is known that for some very low accumulation-rate sites in areas with strong katabatic winds (e.g. “megadunes”, Antarctica) extremely deep CZs can occur, which are potentially able to smooth out even decadal-scale temperature variations (Severinghaus et al., 2010). For this, its deepness would need to be of several dozens of metres, which is highly unrealistic even for glacial Summit conditions (Guillevic et al., 2013; see discussion in Annex A4, p. 1042) as well as for the rather stable Holocene period in Greenland for which no low accumulation and strong katabatic wind situations are to be expected.