How precipitation intermittency sets an optimal sampling distance for temperature reconstructions from Antarctic ice cores

Abstract

Many palaeoclimate proxies share one challenging property: they are not only driven by the climatic variable of interest, e.g. temperature, but they are also influenced by secondary effects which cause, among other things, increased variability, frequently termed noise. Noise in individual proxy records can be reduced by averaging the records, but the effectiveness of this approach depends on the correlation of the noise between the records and therefore on the spatial scales of the noise-generating processes. Here, we review and apply this concept in the context of Antarctic ice-core isotope records to determine which core locations are best suited to reconstruct local- to regional-scale temperatures. Using data from a past-millennium climate model simulation equipped with stable isotope diagnostics we intriguingly find that even for a local temperature reconstruction the optimal sampling strategy is to combine a local ice core with a more distant core $\sim$ 500–1000 km away. A similarly large distance between cores is also optimal for reconstructions that average more than two isotope records. We show that these findings result from the interplay of the two spatial scales of the correlation structures associated with the temperature field and with the noise generated by precipitation intermittency. Our study helps to maximize the usability of existing Antarctic ice cores and to optimally plan future drilling campaigns. It also broadens our knowledge of the processes that shape the isotopic record and their typical correlation scales. Finally, many palaeoclimate reconstruction efforts face the similar challenge of spatially correlated noise, and our presented method could directly assist further studies in also determining optimal sampling strategies for these problems.

1 Introduction

The oxygen and hydrogen isotopic composition of firn and ice recovered from polar ice cores is a key proxy for past near-surface atmospheric temperature changes . Although the physical mechanisms that link local changes in temperature to the isotopic composition of precipitated snow are generally well understood and can be modelled with general circulation models , the quantitative interpretation of ice-core isotope variability, in terms of temperature variability, is complicated by second-order processes that influence the isotopic record, adding noise .

Specifically, the isotopic record that is derived from an ice core is the result of a chain of processes: (1) atmospheric temperature changes along with (2) isotopic fractionation during the pathway from atmospheric moisture to precipitation, (3) the effect of variable and intermittent precipitation, and finally (4) local depositional and post-depositional effects. As we outline in the following, each element of this chain can be associated with a typical spatial length scale over which it is correlated.

Atmospheric temperature variations drive the isotopic composition fractionation of the atmospheric moisture along its pathway to the final stage of precipitation . The spatial coherence of the temperature-related isotopic signal in precipitation is hence determined by the spatial coherence of the variations of the atmospheric temperature field itself. Typical spatial decorrelation scales for temperature anomalies are on the order of $\mathit{\gtrsim }\mathrm{1000}$ km , which implies that ice cores distributed on spatial scales below $\sim \mathrm{1000}$ km should typically record a similar, i.e. correlated, temperature signal. However, the temporal variability of the isotopic composition in the local atmospheric moisture also depends on the variability of the atmospheric circulation, since different air masses may exhibit different source regions and distillation pathways . In addition, the isotopic composition profile across a deposited layer of snow will not directly reflect the temporal variability of the atmospheric isotopic signal due to the intermittent nature of precipitation . By this, the initial isotope signal is weighted with the amount of precipitation, which introduces bias and adds additional variability to the isotopic record . The latter two processes are linked to atmospheric dynamics, and their typical spatial scales range from the mesoscale (i.e. tens of kilometres), driven by topography and orographic effects, to synoptic scales of hundreds of kilometres associated with cyclonic activity and the movement of high- and low-pressure systems. Finally, in polar conditions, the precipitated snow does not directly settle but is constantly eroded, blown away, and redeposited. These depositional processes have been shown to give rise to stratigraphic noise in the isotopic record , which exhibits a small-scale decorrelation scale of a few metres . We further note that the final isotopic record is also influenced by potential exchange processes at the surface and by densification and diffusion within the snow and ice, which are, however, not within the scope of this article.

Both the effect of precipitation intermittency and stratigraphic noise constitute a significant relative contribution to the overall isotopic variability in the form of noise: around a deep drilling site in Dronning Maud Land, East Antarctica, stratigraphic noise was shown to amount to approximately 50 % of the total variance at the seasonal timescale , but quantitative estimates for other Antarctic regions are still missing. A similarly high relative contribution is expected from precipitation intermittency , which probably has a larger impact further inland than compared to coastal regions .

The hierarchy of the different spatial scales of the processes influencing an isotope record determines the effectiveness of reducing the overall noise, since a reduction in the noise level by averaging records will depend on the spatial correlation scale of the different noise sources. For example, if an isotope record were only shaped by temperature variations and stratigraphic noise, it would be sufficient to average records spaced only tens of metres apart, as this would ensure highly correlated temperature signals but uncorrelated stratigraphic noise between the records. However, comparing the correlation-based signal-to-noise ratios derived from nearby isotope records with the signal-to-noise ratios estimated from analysing the records' temporal variability shows that the reproducibility on a local scale does not necessarily imply a climatic, i.e. temperature-driven, origin. Instead, the additional noise sources from circulation variability and precipitation intermittency are likely to exhibit larger decorrelation lengths than the stratigraphic noise . Taking this into account, we expect there to be an optimal length scale which lies between the decorrelation scales of the local noise and of the temperature and which results in a trade-off between averaging out atmospheric circulation and precipitation intermittency effects, while also ensuring a sufficient coherence in the recorded temperature signal.

The aim of the present study is to use data from a climate model equipped with stable isotope diagnostics to systematically study the different typical process scales – including those from atmospheric temperature variations, circulation variability, precipitation intermittency, and the isotope–temperature relationship – to determine the optimal spatial arrangement of ice-core locations which maximizes the correlation with temperature at a specific target site. To address this problem we focus on target sites on the East Antarctic Plateau. Our results show that the average of multiple ice-core isotope records yields a higher degree of correlation with temperature when the sampled locations are spread across distances of 1000 km or more from the target site than when they are all located close ($<\mathrm{250}$ km) to the target site. While these results may seem counterintuitive at first, we qualitatively explain their general features with a simple analytical model that uses the typical spatial correlation structures associated with the temperature and isotope fields, as well as with the noise generated by precipitation intermittency.

2 Data and methods

2.1 Climate model data

We use data from the past-millennium simulation (800–1999 CE; ) of the fully coupled ECHAM5/MPI-OM-wiso atmosphere–ocean general circulation model equipped with stable isotope diagnostics . This simulation is forced by greenhouse gases, volcanic aerosols, total solar irradiance, land use changes, and changes in the Earth's orbital parameters. The model's atmospheric component ECHAM5-wiso is run with a T31 spectral resolution (3.75${}^{\circ }$ $\times$ 3.75${}^{\circ }$) and with $\mathrm{19}$ vertical levels . Compared to observations, the climatological relationship between temperature and the precipitation isotopic composition is reproduced well by the model, but it is biased towards warm temperatures in the T31 setup and its isotopic composition is not depleted enough over Antarctica . These issues can be improved upon by using a higher spatial resolution ; however, such a higher-resolution model is not needed for our study, since we are mainly interested in the relative variability between sites and not in the absolute temperature or isotope values. The full atmosphere–ocean model was compared to observational data and palaeoclimate records for two equilibrium simulations under pre-industrial and Last Glacial Maximum conditions , and the past-millennium simulation was used to reconstruct North Atlantic atmospheric circulation in combination with ice-core isotope data .

In this study, we use the 1200-year ECHAM5/MPI-OM-wiso time series of 2 m surface air temperature ($T_{\mathrm{2}\mathrm{m}}$), precipitation ($p$), and oxygen isotopic composition in precipitation (the relative abundance of oxygen-18 to oxygen-16 isotopes, denoted as ${\mathit{\delta }}^{\mathrm{18}}$O) extracted from the total number of $\mathrm{442}$ model grid cells that are available for the Antarctic continent .

2.2 Data processing

The model simulation output has a monthly temporal resolution, while ice-core isotope records typically exhibit an annual (or even lower) resolution. The latter is commonly achieved by averaging the isotopic data across annual layers of snow and ice, which are determined through a dating approach. The resulting annual isotopic composition data therefore include a weighting effect due to the intra-annual variability in the amount of precipitation. To account for this, we produce two versions of annual data from the monthly model output : (1) the 2 m temperature and oxygen isotopic composition data averaged to an annual resolution without any weighting (denoted as $T_{\mathrm{2}\mathrm{m}}$ and ${\mathit{\delta }}^{\mathrm{18}}$O in the following) and (2) the respective monthly data averaged to an annual resolution including the weighting by the monthly precipitation amount (denoted as precipitation-weighted data $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$).

In extremely dry areas with very little precipitation or high evaporation, numerical instabilities can occur for the modelled isotopic composition in precipitation, resulting in anomalously strong positive or negative spikes in the isotope time series, which is also observed for the Antarctic data in our model simulation. We set a threshold of $\mathrm{4}$ times the interquartile range of a time series, above or below which data points are regarded as outliers, and apply it to every grid cell in order to filter outliers in the ${\mathit{\delta }}^{\mathrm{18}}$O and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series. This approach removes $\mathrm{443}$ anomalous annual values ($<\mathrm{0.1}$ %), out of which $\mathrm{435}$ anomalies occur for the model year $\mathrm{970}$ CE.

2.3 Data analyses

2.3.1 General approach

The overarching aim of this study is to determine a set of locations from which the averaged model data optimally reconstruct the $T_{\mathrm{2}\mathrm{m}}$ temperature time series at a target site, i.e. a specified model grid cell of interest. The optimal reconstruction is assessed by maximizing the Pearson correlation coefficient ($r$) with the target site temperature. To define a spatial set, we combine a given number, $N_{\mathrm{\ell }}$, of model grid cells and vary $N_{\mathrm{\ell }}$ and the distances of these locations relative to the target site. To derive implications for actual ice-core studies, we use the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series at the locations as a surrogate for ice-core isotope records. We thus neglect stratigraphic noise and any further depositional or post-depositional effects on the isotopic record, and therefore our results represent an upper limit of the extent to which ice cores can reconstruct the climatic temperature signal in the atmosphere. In order to learn how the different underlying processes affect the results and to isolate their contributions, we compare the results obtained for ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ with those obtained for $T_{\mathrm{2}\mathrm{m}}$, $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and ${\mathit{\delta }}^{\mathrm{18}}$O. In addition to using only a single target site, we analyse several adjacent target sites in a given region to derive results that are relevant on local to regional spatial scales. In the next section, we present the two main methods that we use to assess the optimal reconstructions.

2.3.2 Assessing optimal reconstructions

Selecting optimal sites

In a first approach, we select an optimal set of ice-core locations to reconstruct a target site's $T_{\mathrm{2}\mathrm{m}}$ time series by sampling without replacement $N_{\mathrm{\ell }}$ grid cells that lie within a selection circle of $\mathrm{2000}$ km radius around the target site and then correlating the average ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series from these $N_{\mathrm{\ell }}$ grid cells with the target site temperature. We perform this for different $N_{\mathrm{\ell }}$ and determine the optimal set of cores for each value of $N_{\mathrm{\ell }}$ from the maximum correlation value across all selection trials. For this, we either sample all possible combinations of grid cell locations within the selection circle, if the number of possibilities does not exceed ${\mathrm{10}}^{\mathrm{7}}$, which effectively applies to all $N_{\mathrm{\ell }}\le \mathrm{3}$, or we randomly sample ${\mathrm{10}}^{\mathrm{7}}$ times from all the possible combinations.

Optimal sampling structure

In order to learn about the typical spatial scales associated with the processes that contribute to the overall temperature–isotope relationship, we aim to investigate how the reconstruction quality depends on the radial distances between the target site and the locations of the ice-core network only, neglecting their angular positions. To do so, one could use the above random selection trials and bin them according to the distances of the selected locations relative to the target site. However, for $N_{\mathrm{\ell }}>\mathrm{1}$ the number of possible grid cell combinations quickly becomes much larger than the actual number of grid cells. In combination with the limited computation time, such an approach would likely result in uneven sample sizes for the available distance combinations for larger $N_{\mathrm{\ell }}$, especially for distances farther away from the target site due to the radially increasing number of grid cells.

Figure 1

Conceptual sketch of the ring sampling approach. Around a given Antarctic target site (black crosses in a and b) we define consecutive rings of $\mathrm{250}$ km radial width (red lines in a and b). From the array of available model grid cells (grey points in b), we choose sets of grid cells which consist of $N_{\mathrm{\ell }}$ cells and which are drawn from $N_{\mathrm{\ell }}$ radial bins determined by a selected combination of rings. As an example for $N_{\mathrm{\ell }}=\mathrm{2}$, possible grid cell sets are shown for the cases of (i) combining the innermost ring with itself (grid cells marked black), (ii) combining the innermost ring with the second ring (grid cells marked blue), and (iii) combining the third and the fourth ring (grid cells marked orange). Also shown in (a) are our main study regions (black polygons) around the EDML (upward-pointing triangle) and Vostok (downward-pointing triangle) ice-core sites. The ring width of $\mathrm{250}$ km is chosen as a trade-off between high spatial resolution and the requirement that a sufficient number of grid cells lie inside each ring. Note that for aesthetic reasons, only four rings are displayed instead of the actually used nine rings and that the model grid is shown simplified as a regular grid in space.

[Figure omitted. See PDF]

Here, we instead use a second more general approach that ensures constant sampling of the entire available space of radial distance combinations and which also reduces local effects in the climate model data and provides more stable correlation results. For a given target site, we define as sampling regions nine concentric rings around the target site with increasing radius in steps of $\mathrm{250}$ km (Fig. ) and identify all grid cells that lie within each of these rings. The sampling of $N_{\mathrm{\ell }}$ grid cells is then implemented in the following two-step process: first, we determine all possible combinations of selecting $N_{\mathrm{\ell }}$ rings with replacement. For every ring combination, we then apply the following second step: we sample one individual grid cell from each of the $N_{\mathrm{\ell }}$ rings (see the examples in Fig. b for an illustration), extract from this grid cell set the time series for a studied model variable, average them, and compute the degree of correlation of this average record with the target site temperature. This second step is iterated over the available number of grid cell sets, and we report the mean correlation across all analysed grid cell sets. For the iteration, we identify all possible grid cell sets until $N_{\mathrm{\ell }}=\mathrm{2}$; for $N_{\mathrm{\ell }}\ge \mathrm{3}$, we resort to Monte Carlo sampling of the grid cell sets due to computational reasons, for which we estimated ${\mathrm{10}}^{\mathrm{5}}$ iterations to provide sufficient convergence of the results.

This approach provides insight into the average spatial structure of the correlation with the target site temperature for sampling $N_{\mathrm{\ell }}$ locations from the model field depending on the radial distances of the locations, as given by the respective ring midpoint radii. We denote this quantity as the sampling correlation structure. Note that in the one-dimensional case ($N_{\mathrm{\ell }}=\mathrm{1}$), the sampling correlation structure is identical to what is often called the spatial correlation structure, i.e. the average correlation as a function of radial distance.

2.3.3 Study regions

To derive sampling correlation structures which are representative on a regional scale, we conduct the above analysis for specific regions by successively using each model grid cell in the region as a target site and then averaging the resulting sampling correlation structures across these target sites.

We make use of this approach for two subregions of the East Antarctic Plateau: the Dronning Maud Land (DML) region in the Atlantic sector of the plateau and the Vostok region in the Indian Ocean sector, both of which include existing deep ice-core drilling sites and large arrays of shallower ice and firn cores. We define the DML region as the area of $\pm \mathrm{17.5}{}^{\circ }$ longitude and $\pm \mathrm{5}{}^{\circ }$ latitude around the European Project for Ice Coring in Antarctica (EPICA) DML site (EDML; 75${}^{\circ }$ S, 0${}^{\circ }$ E; Fig. a), consisting of $\mathrm{26}$ model grid cells. This region encompasses the site of the deep EDML ice core and $>\mathrm{50}$ firn and shallow ice cores . For the Vostok region, we choose an identical latitudinal and longitudinal coverage with respect to the Vostok station (78.47${}^{\circ }$ S, 106.83${}^{\circ }$ E; Fig. a), covering $\mathrm{30}$ model grid cells and encompassing the sites of the deep Vostok and Dome C ice cores, several shallower cores , and the new deep drilling project (“Little Dome C”) where an ice core extending back more than 1 million years is envisaged .

3 Results

3.1 Spatial scale of the temperature anomalies and the local temperature–isotope relationship

First, we assess the extent to which a single ice-core record, i.e. the annual isotope time series of an individual grid cell in the model simulation, is representative of the local- and regional-scale variability of the near-surface atmospheric temperature.

Figure 2

Temperature decorrelation lengths and local temperature–isotope relationship across Antarctica. (a) The temperature decorrelation lengths ($\mathit{\tau }$, in kilometres) for each Antarctic model grid cell estimated by fitting an exponential model to the correlation–distance relationship (see Eq. ) obtained from correlating the local annual near-surface $T_{\mathrm{2}\mathrm{m}}$ time series with the respective temperature time series from all other grid cells. Note that only the continental grid cells were used for the fit. Although the decorrelation lengths show a strong partition between East and West Antarctica, they are larger than $\mathrm{1000}$ km at most locations. (b, c) The local correlations at each model grid cell between the annual time series of precipitation-weighted oxygen isotope composition and of (b) near-surface temperature and (c) precipitation-weighted near-surface temperature. The difference between the maps (d) clearly demonstrates that precipitation intermittency is a major limiting factor for the temperature–isotope relationship.

[Figure omitted. See PDF]

The temperature field over Antarctica in the climate model exhibits large-scale coherent variations (Fig. a) with a clear two-part structure, which is roughly divided by the Transantarctic Mountains: for most parts of the East Antarctic Plateau, the temperature field shows typical decorrelation lengths between $\sim \mathrm{1500}$ and $\mathrm{2500}$ km, while the decorrelation lengths are notably lower, with values ranging from $\sim \mathrm{500}$ to $\mathrm{1500}$ km, for larger parts of the West Antarctic Ice Sheet and for the Antarctic Peninsula. Still, for perfect ice cores, i.e. assuming an ideal temperature proxy record that is only governed by local temperature variations, a single ice core would capture the temperature variability in both East and West Antarctic regions across hundreds of kilometres.

However, as simulated by the isotope-enabled climate model, actual single Antarctic ice-core isotope records only explain a low portion of the variations in the local temperature fields: correlating the annual precipitation-weighted field of modelled ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ with the annual $T_{\mathrm{2}\mathrm{m}}$ time series at the same grid cell results in generally low correlations (mean of $\mathrm{0.38}$), which across all analysed grid cells range from $\sim \mathrm{0.1}$ up to $\sim \mathrm{0.57}$, with $\sim \mathrm{60}$ % of the correlations $\le \mathrm{0.4}$ (Fig. b). The correlations are overall improved when the $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ time series is used instead of the $T_{\mathrm{2}\mathrm{m}}$ time series (mean correlation of $\mathrm{0.53}$, range $\sim \mathrm{0.1}$ to $\mathrm{0.77}$; Fig. c) but with unaffected correlation values mostly in the coastal regions (Fig. d). This shows that precipitation intermittency is a major limiting factor for the local temperature–isotope correlation on the continental plateau but is less important near the coasts due to higher and more regular snowfall amounts there .

3.2 Spatial correlation with local temperature

In the next step, we investigate how a local temperature record correlates in space with the temperature itself and with the oxygen isotope composition. For this, we choose the EDML and Vostok drilling sites as target sites and correlate the annual $T_{\mathrm{2}\mathrm{m}}$ time series at these target sites with the spatial fields of annual temperature and of annual ${\mathit{\delta }}^{\mathrm{18}}$O, both unweighted and weighted by the precipitation amount (Fig. ).

Figure 3

Spatial correlation with the temperature at the EDML and Vostok target sites. Shown are the correlations of the $T_{\mathrm{2}\mathrm{m}}$ time series at the target sites EDML (a–d) and Vostok (e–h) with the spatial fields of temperature (a, d), precipitation-weighted temperature (b, f), oxygen isotope composition (c, g), and precipitation-weighted oxygen isotope composition (d, h). The target sites are marked with a black cross; black lines indicate correlation contour lines incremented in steps of 0.1.

[Figure omitted. See PDF]

We find that the correlation patterns with the temperature field itself are largely radially symmetric with respect to the target sites and decay uniformly with distance within the first couple of hundred kilometres from the target (Fig. a, e). However, for ${\mathit{\delta }}^{\mathrm{18}}$O, and also partly through the effect of the precipitation weighting, radial asymmetry in the correlation patterns occurs. This is particularly striking for the EDML target site. Here, the maximum in correlation with the ${\mathit{\delta }}^{\mathrm{18}}$O field is not centred on the target site but displaced by $\sim \mathrm{1200}$ km towards the southeast (Fig. c, d). Some spatial displacement in maximum correlation is also visible for the Vostok target site and the $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, ${\mathit{\delta }}^{\mathrm{18}}$O, and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ fields (Fig. f–h), but in different directions between $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ and the oxygen isotope fields and much smaller than in the case of EDML. We also note that the correlation patterns for the $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, ${\mathit{\delta }}^{\mathrm{18}}$O, and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ fields still contain radially symmetric contributions with respect to the target sites, which are more pronounced for the Vostok than for the EDML target site.

3.3 Selecting optimal ice-core sites for temperature reconstructions

The above analyses have shown firstly that isotope records from single ice cores likely only capture a small portion of the local interannual temperature variability, suggesting that additional processes, such as precipitation intermittency, influence the isotopic signal and decrease the degree of correlation with the local temperature record. Interpreting these additional processes as noise raises the question of whether the correlation with temperature can be improved upon by averaging isotope records across space. In addition, we have seen that the correlation of an oxygen isotope composition record with a local temperature record is not necessarily maximal at the location of the temperature recording, posing the question of how locations of isotope records should be spatially arranged with respect to the location of the temperature record in order to get the best correlation. To address these questions, we assume an ideal world in which the climate model data are a perfect surrogate for the true climate and proxy variations at each site, and we set up the simple experiment of selecting and averaging ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ records from grid cells within a $\mathrm{2000}$ km circle around a target site (see Sect.  for details) to determine what spatial array of $N_{\mathrm{\ell }}$ ice cores optimizes the temperature correlation with the target site.

Figure 4

Selecting ice-core locations that optimally reconstruct interannual temperatures at the EDML and Vostok drilling sites. The maps show the correlation coefficient in the climate model data between the annual temperature time series at the target sites (black crosses) EDML (a–c) and Vostok (b–f) with the time series fields of precipitation-weighted oxygen isotope composition (${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$). Filled black points denote grid cells that yield the maximum correlation between the target site temperature and the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series from either selecting a single grid cell ($N_{\mathrm{\ell }}=\mathrm{1}$; a, d) or from averaging across $N_{\mathrm{\ell }}=\mathrm{3}$ (b, e) or $N_{\mathrm{\ell }}=\mathrm{5}$ (c, f) grid cells, obtained from iteratively selecting sets of $N_{\mathrm{\ell }}$ grid cells from within a selection circle of $\mathrm{2000}$ km radius around the target site indicated by the black radial lines (see Sect.  for details). Interestingly, non-local ice-core locations systematically show the strongest relationship with the target site temperature.

[Figure omitted. See PDF]

For our specific model simulation and specifying the EDML drilling site as the target site, we already know from Fig. a that the optimal location for a single ice core is not the local grid cell at the target site but should be a $\sim \mathrm{1200}$ km southeastward site. Choosing this more distant site increases the correlation with the target temperature from an $r$ value of $\mathrm{0.30}$ for the local EDML site to a value of $\mathrm{0.44}$ (Fig. a). Even more intriguingly, when we analyse the maximum correlations with the EDML target temperature for an average of three or five cores chosen from the $\mathrm{2000}$ km selection circle (Fig. b–c), we find optimal locations that in both cases are scattered at significant distances around the target and which yield an even further increase in correlation ($r=\mathrm{0.50}$ for $N_{\mathrm{\ell }}=\mathrm{3}$, $r=\mathrm{0.52}$ for $N_{\mathrm{\ell }}=\mathrm{5}$). We obtain comparable results when the Vostok drilling site is specified as the target (Fig. d–f). The optimal single core would be at a location $\sim \mathrm{190}$ km west of Vostok ($r=\mathrm{0.49}$ compared to the local correlation of $r=\mathrm{0.46}$), and the optimal locations for averaging three or five cores again all lie scattered around the target without including it and result in a significant further increase in correlation ($r=\mathrm{0.60}$ for $N_{\mathrm{\ell }}=\mathrm{3}$, $r=\mathrm{0.63}$ for $N_{\mathrm{\ell }}=\mathrm{5}$).

We generalize these findings by considering each Antarctic model grid cell as a target site and determining in each case the ice-core location that results in an optimal correlation with the target site. As in the above EDML case study, about half of the optimal locations for a single ice core are situated at distances between $\mathrm{500}$ and $\mathrm{1500}$ km from the respective target sites, while only about $\mathrm{10}$ % lie within $\mathrm{500}$ km from the targets. We note that this distribution might be affected by the number of available sampling points (i.e. model grid cells) per distance bin, which increase with increasing distance from the target site. However, after weighting the distance distribution with the average inverse number of available grid cells per distance bin, still only about one-fifth of the optimal distances lie within $\mathrm{500}$ km from the targets.

3.4 Optimal ice-core sampling structures

The approach for choosing optimal ice-core locations yields straightforward and instructive results. However, it might be doubtful as to whether these results can be directly applied to the real world, since they might depend on the specific simulated climate state, depend on the specific climate model and model isotope scheme used, or result from statistical overfitting. We therefore adapt our approach in a next step to learn more about the general spatial arrangement of the optimal ice-core locations which yield the maximum correlation with temperature. This is done by applying our concept of sampling correlation structures (see Sect.  and the illustration in Fig. ), which studies the correlation patterns only as a function of radial distance from the target site by averaging across $\mathrm{250}$ km radial bins and across the angular positions, thereby reducing local variability in the model data. Additionally, we apply the approach to all target sites in our DML and Vostok study regions (Sect. ) and average the results across these sites to obtain regional estimates. Finally, we analyse each of the model variables to highlight the differences between the individual fields.

Figure 5

Average sampling correlation structures with temperature for the DML and Vostok regions in the case of sampling single locations. Shown as a function of distance is the average correlation between the interannual near-surface temperature ($T_{\mathrm{2}\mathrm{m}}$) at a target site and the spatial fields of $T_{\mathrm{2}\mathrm{m}}$ (black), oxygen isotope composition (${\mathit{\delta }}^{\mathrm{18}}$O, green), and precipitation-weighted oxygen isotope composition (${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$, blue). Averaging was performed in two steps: first, for a given target site, the correlations with the target site temperature were averaged across grid cells lying within $\mathrm{250}$ km wide consecutive rings around the given target site. Secondly, this analysis was conducted for all target sites in the DML (a) and Vostok (b) region, and the results were averaged across the respective region (see Sect.  and for details). Shading denotes $\pm \mathrm{1}$ standard deviations of the correlation results across the different target sites in each region. The black dashed lines indicate an exponential fit to the $T_{\mathrm{2}\mathrm{m}}$ data.

[Figure omitted. See PDF]

When we sample only a single location ($N_{\mathrm{\ell }}=\mathrm{1}$), the sampling correlation structure is conceptually equivalent to the average correlation with distance, and it therefore simply gives the spatial decorrelation in the case of sampling from the $T_{\mathrm{2}\mathrm{m}}$ field itself. The average sampling correlation structures for $T_{\mathrm{2}\mathrm{m}}$ across the DML and Vostok regions (Fig. ) can be described by an exponential decay with a length scale of $\sim \mathrm{1900}$ km in both cases, consistent with the estimated spatial temperature decorrelation lengths for the individual grid cells in these regions (Fig. a). In accordance with the general expectation, the maximum average correlation with the target site temperature is thus obtained from sampling the innermost ring only.

When we compare these results to the average sampling correlation structure for the ${\mathit{\delta }}^{\mathrm{18}}$O field, we find for the DML region a much lower average correlation with the target site temperature as a function of distance (Fig. a). The average correlation for the innermost ring ($<\mathrm{250}$ km) is $\sim \mathrm{0.4}$, but it decreases only slightly within the first $\sim \mathrm{800}$ km, followed by a slightly steeper decrease and nearly constant levels of $r\mathit{\lesssim }\mathrm{0.2}$ for distances $\mathit{\gtrsim }\mathrm{1600}$ km. For the Vostok region (Fig. b), the average sampling correlation structure for ${\mathit{\delta }}^{\mathrm{18}}$O exhibits a nearly linear decrease from an initial value of $r\sim \mathrm{0.6}$ to $r\sim \mathrm{0.1}$ in the final ring ($>\mathrm{2000}$ km). When we analyse the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ fields we find that precipitation weighting reduces the correlation values in both regions but that it does not have a large effect on the shape of the sampling correlation structures itself.

Extending this analysis to the two-dimensional case of sampling and averaging $N_{\mathrm{\ell }}=\mathrm{2}$ locations offers the possibility of investigating the average correlation not only as a function of distance from the target site but also implicitly as a function of distance between the two sampled locations (Fig. ). The difference in the average sampling correlation structure between the fields of $T_{\mathrm{2}\mathrm{m}}$ and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ is even more pronounced for $N_{\mathrm{\ell }}=\mathrm{2}$ than for $N_{\mathrm{\ell }}=\mathrm{1}$. The maximum average correlation for $T_{\mathrm{2}\mathrm{m}}$ is still found when both sampling locations lie inside the innermost ring, as shown for the DML region (Fig. a). However, for ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ the optimal arrangement of two locations is to sample one location from within the innermost ring but the second location from within the fifth ring, i.e. between $\sim \mathrm{1000}$ and $\mathrm{1250}$ km from the target site (Fig. c). Part of this structure is related to the effect of precipitation intermittency, which can be seen from the average sampling correlation structure of the $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ field (Fig. b). Here, in contrast to $T_{\mathrm{2}\mathrm{m}}$, the correlation is about as high when we combine the innermost ring and one ring further away as when we sample both locations from within the innermost ring.

Figure 6

Average sampling correlation structure with temperature for the DML and Vostok regions in the two-dimensional case of sampling two locations. Shown is the mean correlation of all possible single correlations between the target site temperature and the average of two grid cells of (a, d) $T_{\mathrm{2}\mathrm{m}}$, (b, e) $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and (c, f) ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series sampled from the same ring or from two different rings. This analysis was conducted for every target site in the DML region (a–c) and in the Vostok region (d–f), and the results were then averaged across the respective region. For each plot, the axes display the distance from the target site; the $x$ ($y$) axis represents the first (second) sampled ring, with the results being mirrored along the diagonal for aesthetic reasons. The tick marks indicate the border distances of the rings. Note the marked difference in the locations of the correlation maxima between $T_{\mathrm{2}\mathrm{m}}$ and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ for the DML region, and also for the Vostok region the – albeit marginal – correlation maximum for ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ is achieved by combining the innermost ring with the ring between $\mathrm{500}$ and $\mathrm{750}$ km.

[Figure omitted. See PDF]

Analysing the Vostok study region leads to comparable results (Fig. d–f), with a similar difference in average sampling correlation structure between $T_{\mathrm{2}\mathrm{m}}$ and $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ as for the DML region and a similar structure of $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ for distances $\mathit{\lesssim }\mathrm{1000}$ km. However, the results for the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ field (Fig. f) do not display such a pronounced maximum correlation when one location is sampled from within the innermost ring and the second one from inside a ring further away as is observed for the DML region. This suggests that the regional differences in the spatial correlation structure of the ${\mathit{\delta }}^{\mathrm{18}}$O field (Fig. ) have an influence here.

The general feature of the optimal ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ sampling arrangement is robust throughout Antarctica, despite the above regional differences. When we analyse all available Antarctic target sites and fix the first core location to lie inside the innermost ring, in $\sim \mathrm{82}\mathit{\%}$ of all cases the optimal second core location is at least the second ring ($>\mathrm{250}$ km), and in $\sim \mathrm{63}$ % of the cases it is the second to fourth ring ($\mathrm{250}$–$\mathrm{1000}$ km).

Figure 7

The optimal arrangement for averaging three or five ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ ice cores to reconstruct the target site temperature at the EDML (a, c) and Vostok (b, d) drilling sites. Displayed are subsets of the sampling correlation structures for $N_{\mathrm{\ell }}=\mathrm{3}$ and $\mathrm{5}$, showing along the vertical axis the optimal five of all possible combinations of rings (best denoted as rank 1, fifth best as rank 5), i.e. those which exhibit the five highest mean correlation values across ${\mathrm{10}}^{\mathrm{5}}$ random trials of averaging $N_{\mathrm{\ell }}=\mathrm{3}$ (a, b) or $N_{\mathrm{\ell }}=\mathrm{5}$ (c, d) grid cells from these rings. The ring bin borders are marked by thin vertical lines with their distances from the target site given on the horizontal axes; the selected optimal ring combinations are marked as black dots. Systematically, arrangements which combine ice cores from the innermost ring with ice cores further away are found to be optimal, with larger distances for the EDML target site.

[Figure omitted. See PDF]

We also obtain similar results when averaging $N_{\mathrm{\ell }}=\mathrm{3}$ or $\mathrm{5}$ locations of the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ field to reconstruct the target site temperature (Fig. ). For computational reasons, we only analyse single target sites here. When EDML is set as the target site, the optimal sampling configuration is such that one to two core locations lie in the innermost ring, while the others are distributed at distances mostly between $\sim \mathrm{750}$ and $\mathrm{1500}$ km from the target. For reconstructing the Vostok target site temperature, the optimal core locations combine the innermost ring with locations distributed mostly across the second to third ($\mathrm{250}$–$\mathrm{750}$ km) ring.

Figure 8

Gain in correlation and risk of adverse sampling. (a) The average correlation with the target temperature at the EDML (red) and Vostok (blue) sites depending on the number of locations, $N_{\mathrm{\ell }}$, used for averaging the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series. Sampling is performed either locally from the innermost ring only (dashed lines) or from all possible individual combinations of locations for the respective optimal ring combination determined for each $N_{\mathrm{\ell }}$ (solid lines). Compared to the local samples, which show virtually no or only a small increase with the number of sampled locations, the correlation increases markedly with $N_{\mathrm{\ell }}$ when sampling from the optimal rings, as highlighted by the shaded area. (b) Histogram of individual correlations for sampling from the optimal ring combination when averaging $N_{\mathrm{\ell }}=\mathrm{3}$ locations compared to the correlation (vertical lines) for sampling from the innermost ring only, displayed for the EDML (red) and Vostok (blue) target sites. In both cases, the correlation is higher than the local value for more than 93 % of the optimal ring combination samples.

[Figure omitted. See PDF]

In summary, averaging the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ time series across the optimal locations clearly increases the average correlation with the target site temperature more strongly with the number of locations compared to sampling all core locations only locally close to the target site, i.e. from the grid cells that lie within the innermost ring (Fig. a). While the local correlation for the EDML target site stays constant around $\mathrm{0.31}$, the optimal correlation rises to $\mathrm{0.35}$ for $N_{\mathrm{\ell }}=\mathrm{2}$ and to $\mathrm{0.43}$ for $N_{\mathrm{\ell }}=\mathrm{10}$, which is equivalent to nearly a doubling in the explained variance. For the Vostok target site, we observe a nearly concurrent increase in correlation between the local and optimal sampling up until $N_{\mathrm{\ell }}=\mathrm{2}$ from $\mathrm{0.45}$ to $\sim \mathrm{0.50}$, but for larger $N_{\mathrm{\ell }}$ the optimal correlation also increases more strongly and reaches $\mathrm{0.58}$ for $N_{\mathrm{\ell }}=\mathrm{10}$, which is a $\sim \mathrm{1.7}$-fold higher explained variance compared to $N_{\mathrm{\ell }}=\mathrm{1}$.

These results are the mean value from averaging across many possible combinations of individual locations. In reality, any new drilling campaign or reanalysis of existing ice cores only represents one single combination of locations. Therefore, we further assess the risk of an “adverse optimal sampling”, i.e. the probability of choosing by chance a specific sampling realization from the optimal ring combination which yields a lower correlation than the correlation for sampling locally. For this purpose, we compare the distribution of individual correlations from sampling the optimal ring combination with the value obtained from sampling only the local sites which lie in the innermost ring. Overall we find the risk of adverse optimal sampling to be low, since more than 93 % of all individual correlation values in the example of $N_{\mathrm{\ell }}=\mathrm{3}$ are actually larger than the respective local correlation (Fig. b).

4 Discussion

4.1 Dependence on radial distance

Oxygen isotope records derived from ice cores are commonly interpreted to reflect local temperature changes at the ice-core drilling site. Here we have shown that while there is local isotope–temperature correlation (Fig. b), this correlation can be increased considerably by averaging isotope records across space (Fig. a) following a distinct radial pattern which combines the local target site with locations between a few hundred kilometres and $\sim \mathrm{1000}$ km from the target site (Figs. c, f, and ). These results are based on a method which investigates the spatial correlation structure only as a function of radial distance by averaging across the azimuthal component. The motivation for this approach is that from physical arguments we expect the first-order spatial correlation patterns to be invariant against rotation. Such radial symmetry is indeed observed as the leading component of the spatial correlation structure of the temperature field and as a second-order component of the oxygen isotope field (Fig. ). We interpret these symmetric contributions as a general feature of the underlying atmospheric processes compared to individual, local correlation maxima which are more due to the actual dynamics. Therefore, we expect that our results obtained from the radial sampling correlation structures should be largely independent of the climate state, or the climate model used, and thus serve as valid recommendations for real-world applications. In the next section, we substantiate this interpretation by showing that a simple conceptual model can predict the sampling correlation structure from the basic processes which shape the isotopic composition time series modelled only as a function of radial distance. Finally, we will discuss the relevance of our results to actual ice-core studies.

4.2 Conceptual model of the optimal sampling structure

For a conceptual model of the sampling correlation structure, we focus on the three main atmospheric processes that influence the oxygen isotope records in ice cores: (i) temperature variations, (ii) precipitation intermittency, and (iii) the temperature–isotope relationship. We statistically model the associated fields of $T_{\mathrm{2}\mathrm{m}}$, $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ separately in order to understand the influence of each process (see Appendix  for details), and we assess, for comparable results, the predicted average sampling correlation structure with the target site temperature in the two-dimensional case of averaging two locations in the same manner that we analysed the climate model data.

To model the atmospheric temperature field, we assume an isotropic exponential decay of the spatial correlation with a constant decorrelation length (Appendix ). Such an exponential temperature decorrelation is a commonly observed feature and also confirmed by our climate model data (Figs. a, a, e, and ). Given this relationship, we find good agreement for the two-dimensional sampling correlation structure between the conceptual model and the climate model data regarding both absolute correlation values and the spatial pattern (Fig. a). We emphasize that the maximum correlation with the target site temperature naturally occurs in the case of an isotropic correlation decay when the averaged two (or $N_{\mathrm{\ell }}$) locations are close to the target site, as any location which is further away will contribute a temperature signal that is less similar to the other locations.

To elucidate the role of precipitation intermittency, we follow the simplest assumption, which is that this process can be described by partly aliasing the original temperature signal into temporal white noise . We further assume that this noise is not independent between sites but that it follows the spatial scale of precipitation events, which we describe as an exponential decorrelation in space with a second length scale (Appendix ). This intermittency length scale is related to the atmospheric processes that deliver precipitation, e.g. synoptic systems, and is hence assumed to be smaller than the length scale of the temperature anomalies. The introduction of this second length scale into our conceptual model generally explains the optimal sampling structure we obtained from the climate model data. Qualitatively, close locations exhibit a strong correlation in temperature but also in the noise from precipitation intermittency; therefore, this noise cannot be reduced by averaging the locations, yielding an overall low signal-to-noise ratio. However, with increasing distance between the locations, the intermittency noise decorrelates faster than the temperature field due to the different decorrelation scales, resulting in an optimal distance of maximum signal-to-noise ratio. This is also reflected in our conceptual model (Figs.  and b, e): when fixing the position of one core to the innermost location and varying only the distance from the target site of the second core location, the correlation with the target site temperature first increases with increasing distance of the second location and then maximizes at an optimal distance before it decays with a further increase in distance. In the climate model data, we observed a similar feature for the precipitation-weighted temperature (Fig. ), though it was not as clear as in the conceptual model. This mismatch could be related to the assumed isotropy in the conceptual model and the according azimuthal averaging done in the climate model data analysis, which potentially smears the intermittency effect in the climate model data due to slight differences in the decorrelation lengths between the different horizontal directions.

In order to incorporate the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ field into the conceptual model, we need to account for the spatial temperature–isotope relationship. To accomplish this, we parameterize the spatial dependence of the correlation between temperature and the oxygen isotope composition with a simple isotropic linear model based on the climate model data results (Fig.  and Appendix ). In addition, we assume that the same effect of precipitation intermittency that we adopted for the temperature field is also applicable to the oxygen isotope field. With these simple assumptions, we obtain good qualitative agreement for the DML region between the conceptual model and the climate model data results (see Figs. c and c). In addition, when we change the parameterized isotope–temperature relationship such that it more closely resembles the Vostok region data (Fig. b), the sampling correlation structure in the conceptual model (Fig. f) is more similar to the observed correlation structure (Fig. f). However, in general the conceptual model fails for ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ to reproduce the actual range in correlations as it produces much lower values than expected.

In summary, our conceptual model provides a quantitative understanding of the spatial correlation of the temperature in the climate model data and at least a qualitative understanding of the processes that affect the correlation between temperature and the ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ field, i.e. precipitation intermittency and the spatial temperature–isotope relationship. The deficiencies in the conceptual model may be attributed to its simplicity. For the governing processes, we assumed spatially constant and isotropic length scales, neglecting local and direction-related differences in e.g. temperature decorrelation lengths (see Fig. a) or the spatial extent of the coherence of precipitation intermittency. Instead of being constant, the latter may differ depending on the type of precipitation, e.g. synoptic versus stratiform precipitation, and may exhibit directional dependencies related to topography. Furthermore, we assumed constant variance of all time series, thereby ignoring potential weighting effects on the correlations for the spatial average of several locations due to different variabilities between them.

4.3 Relevance for ice-core studies

Our results from analysing the climate model data provide guidance on where to drill or from where to analyse $N_{\mathrm{\ell }}=\mathrm{1}$, 2, 3, or more ice cores in order to optimally reconstruct the atmospheric temperature signal for a certain target site. For this, our analysis highlights two distinct approaches.

The first possibility is to follow the recommendations obtained from directly choosing the specific locations which maximize the correlation with the target site temperature (Fig. ). This is straightforward; however, for applications such locations would need to be derived for every target site in question. In addition, as outlined above, it is unclear whether the results can be one-to-one transferred to the real world, since they might be due to unaccounted model deficiencies or depend on dynamical processes in the atmosphere which could differ between climate states or depend on initial conditions. One indication for this is that we obtain different optimal single core locations for more than half of all investigated Antarctic target sites when we analyse only the first or only the second half of the respective climate model time series compared to the full 1200 years.

We have argued above that the sampling correlation structures, obtained from averaging the individual correlations across space for combinations of concentric rings around the target site, are a more general quantity, and we have shown with our conceptual model that they are on average governed by the interplay of the different underlying correlation length scales. We expect the latter to vary less between different climate periods or states or between regions. This is substantiated by the fact that the sampling correlation structures for two cores (Fig. ) are much more robust against analysing only the first or the second half of the model time series, which is different to the results from directly choosing optimal locations. Thus, the sampling correlation structure offers a general approach for finding an optimal ice-core network, but with the downside that it informs us only about the relative radial distances of the optimal network around the target site.

Using the sampling correlation structures we arrive at the following recommendations for optimal ice-core sampling configurations. If it is only possible to drill or analyse a single ice core, our results show that it is always best to sample locally, i.e. to place this core near the target site of interest. This is also common practice, given the usual interpretation of ice-core isotope records as a proxy for local temperatures. However, due to the effect of precipitation intermittency modulated by the spatial coherence of the temperature–isotope relationship, it is no longer optimal in the case of drilling two ice cores to collect both cores near the target site, but instead to drill one core at the target site and one at least $\mathrm{500}$ km away. Where three or more ice cores will be drilled or analysed, we expect the optimal spatial configuration to be more dependent on the study region, but our results indicate that it is still likely best to place one core near the target site and distribute the others across several hundred kilometres.

These inferences are based on data from a single climate model simulation together with a simple statistical conceptual model, which should be tested against observations. As a proof of concept, we thus need to create an isotope record that is in first order only governed by temperature variations and precipitation intermittency and remove the impact of small-scale stratigraphic noise from the actual measured records (assuming that any further processes in the pre-depositional to depositional phase contribute negligibly to the local isotopic variations). To accomplish this, one possible strategy would be to use trench sampling campaigns (see , for the EDML site). Then, one test of our optimal sampling configurations could be to combine one trench record, e.g. one from EDML, with another trench sampled at the optimal distance based on our results for $N_{\mathrm{\ell }}=\mathrm{2}$ and correlate the average of these two trench records with the instrumental temperature data set available for EDML. Based on the results in this study we would expect a higher degree of correlation in this case compared to using only one local trench record from EDML. We acknowledge that such an approach would be challenging due to the small amount of available instrumental data ($\sim \mathrm{20}$ years for EDML) and the inevitable dating uncertainties between the two trench records.

Finally, we note that our implications concerning optimal ice-core sampling configurations might in reality be affected by two further processes we have neglected here. Firstly, clear-sky precipitation (“diamond dust”) is a common phenomenon in Antarctica, especially in the drier regions of the Antarctic Plateau, which potentially occurs more regularly than convective-type or stratiform precipitation. Diamond dust formation is not explicitly simulated by the ECHAM5 model, so it is possible that the precipitation intermittency modelled in our simulation is partly offset in reality by a stronger relative contribution of diamond dust to the total precipitation amount. Secondly, surface–atmosphere vapour exchange between precipitation events might constitute a second process which imprints an atmospheric temperature signal into the surface snow, next to precipitation e.g.. This process could hence also partly counteract the impact of precipitation intermittency, depending on its relative importance for the isotopic composition of the surface snow. However, there is no clear consensus in the recent literature on this question, and ultimately we need quantitative estimates of the importance of vapour exchange processes across temporal scales. In any case, these considerations do not affect our general notion that the optimal ice-core sampling configuration depends on the differences in spatial decorrelation scales of the processes which shape the isotopic records.

5 Conclusions

In this study we assessed the spatial sampling configuration of ice cores to optimally reconstruct the annual near-surface temperature at a specific target site. This problem was motivated by the expectation that the major processes influencing the isotopic records of ice cores operate on different spatial scales.

Indeed, by analysing the temperature and isotope data of an isotope-enabled atmosphere–ocean climate model simulating the climatic history over the last millennium in Antarctica, we showed that while in the optimal setup a single ice core should be placed close to the target site of interest, a second core should be located far ($>\mathrm{500}$ km) from the first core. While this may seem surprising at first glance, it can be straightforwardly explained by the interplay of two different correlation lengths in space: one for the temperature anomalies and one parameterizing the spatial coherence of the effect of precipitation intermittency, as demonstrated by a simple conceptual model. Despite the fact that these results were specifically obtained for two regions of the East Antarctic Plateau, we expect similar results to hold for other parts of Antarctica and potentially also for other large-scale ice-coring regions such as Greenland, as long as our simplified assumptions of nearly isotropic exponential decorrelation length scales are also valid there. However, we also suggest verifying our results with a different isotope-enabled climate model in order to rule out any dependence on the specific atmospheric model and isotope scheme applied in the simulation used here.

Overall, our study explicitly improves the planning of drilling or analysis campaigns for spatial networks of ice-core isotope records. In addition, it provides a strategy to analyse an optimal configuration of sampling locations for any proxy which is influenced by two or more processes that exhibit different spatial correlation scales. This likely applies to various marine and terrestrial proxy types, and our strategy might thus offer a step forward in the best use of sampling and measurement capacity for quantitative climate reconstructions, which needs to be investigated in further studies.

Appendix A Conceptual model of sampling correlation structures

A1 General model

We set up a conceptual model for the correlation between a target temperature time series and a spatial average based on a set of locations sampled from a climatic field (sampling correlation structure). Our model assumes simple isotropic and exponential decorrelation structures for the involved climatic fields and is based on previous work which suggests that precipitation intermittency can be described by partly aliasing the original temperature signal into white noise .

In the model, we consider a temperature time series $T_{\mathrm{0}}$ at some target site ${\mathit{r}}_{\mathrm{0}}$ and a scalar field $x$ of a given climate variable. From this field, we select $N_{\mathrm{\ell }}$ time series $x_i$ at the locations ${\mathit{r}}_i$, $i=\mathrm{1},\mathrm{\dots },N_{\mathrm{\ell }}$ and denote the spatial average of these time series by $\stackrel{\mathrm{‾}}{x}=\frac{\mathrm{1}}{N_{\mathrm{\ell }}}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}}{x}_i$. The distances of the $N_{\mathrm{\ell }}$ locations from the target site and the distances between the locations are given by $r_i=|{\mathit{r}}_i-{\mathit{r}}_{\mathrm{0}}|$ and by $d_{ij}=|{\mathit{r}}_i-{\mathit{r}}_j|$, respectively. The correlation between $T_{\mathrm{0}}$ and $\stackrel{\mathrm{‾}}{x}$ follows from

A1$$\mathrm{cor}\left(T_{\mathrm{0}},\stackrel{\mathrm{‾}}{x}\right)={\displaystyle \frac{\mathrm{cov}\left(T_{\mathrm{0}},\stackrel{\mathrm{‾}}{x}\right)}{\sqrt{\mathrm{var}\left(T_{\mathrm{0}}\right)\mathrm{var}\left(\stackrel{\mathrm{‾}}{x}\right)}}},$$ and it is governed by the covariance between the temperature at the target site and the climate field at the sampling locations ${\mathit{r}}_i$, A2$$\mathrm{cov}\left(T_{\mathrm{0}},\stackrel{\mathrm{‾}}{x}\right)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{\ell }}}}\sum _i^{N_{\mathrm{\ell }}}\mathrm{cov}\left(T_{\mathrm{0}},x_i\right),$$ as well as by the covariance between the sampling locations through the variance of their spatial average, A3$$\mathrm{var}\left(\stackrel{\mathrm{‾}}{x}\right)={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{\ell }}^{\mathrm{2}}}}\left(\sum _i^{N_{\mathrm{\ell }}}\mathrm{var}\left(x_i\right)+\mathrm{2}\sum _i^{N_{\mathrm{\ell }}-\mathrm{1}}\sum _j^{N_{\mathrm{\ell }}}\mathrm{cov}\left(x_i,x_j\right)\right).$$ In our model, these quantities depend on the distance between sites and on the correlation structure of the respective field $x$, as we show in the following and as illustrated in Fig. .

Figure A1

Illustration of the decorrelation lengths in the conceptual model. Shown as a function of distance are the correlation between two temperature time series (black), between the intermittency noise (purple), between a temperature and a precipitation-weighted temperature time series (dashed black), between two precipitation-weighted temperature time series (green), and between a target temperature time series and the average of two precipitation-weighted temperature time series (orange) when one is located at the target site and the other one is located away from the target site at a distance as indicated on the $x$ axis. Model parameters are taken from the DML region. The decorrelation curve of the precipitation-weighted temperature time series is simply the superposition of the temperature decorrelation and the decorrelation of the intermittency noise, depending on the intermittency factor $\mathit{\xi }$.

[Figure omitted. See PDF]

A2 Temperature

For the near-surface temperature field, $x\equiv T$, we assume a spatially constant variance, $\mathrm{var}\left(T_{\mathrm{0}}\right)=\mathrm{var}\left(T_i\right)\equiv {\mathit{\sigma }}_T^{\mathrm{2}}$, and an isotropic decorrelation following an exponential decay with a decorrelation length $\mathit{\tau }$; i.e. the covariance between sites is (see black line in Fig. ) $$\begin{array}{}\text{A4}& {\displaystyle }\mathrm{cov}\left(T_{\mathrm{0}},T_i\right)={\mathit{\sigma }}_T^{\mathrm{2}}\exp \left(-{\displaystyle \frac{r_i}{\mathit{\tau }}}\right),\text{A5}& {\displaystyle }\mathrm{cov}\left(T_i,T_j\right)={\mathit{\sigma }}_T^{\mathrm{2}}\exp \left(-{\displaystyle \frac{d_{ij}}{\mathit{\tau }}}\right).\end{array}$$ The correlation between the target site temperature and the spatial average of $N_{\mathrm{\ell }}$ temperature time series is then obtained from

A6$$\mathrm{cor}\left(T_{\mathrm{0}},\stackrel{\mathrm{‾}}{T}\right)={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}}\exp \left(-\frac{r_i}{\mathit{\tau }}\right)}{\sqrt{N_{\mathrm{\ell }}+\mathrm{2}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}-\mathrm{1}}{\sum }_{j=i+\mathrm{1}}^{N_{\mathrm{\ell }}}\exp \left(-\frac{d_{ij}}{\mathit{\tau }}\right)}}}.$$

A3 Precipitation-weighted temperature

To model the effect of precipitation intermittency, we follow and assume that precipitation intermittency redistributes the energy of the temperature time series constantly across frequencies, i.e. creating temporal white noise without changing the total variance. Then, the precipitation-weighted temperature time series at location ${\mathit{r}}_i$ arises from $T_i$ as

A7$$T_i^{\left(\mathrm{pw}\right)}={\left(\mathrm{1}-\mathit{\xi }\right)}^{\mathrm{1}/\mathrm{2}}{T}_i+{\mathit{\xi }}^{\mathrm{1}/\mathrm{2}}{\mathit{\sigma }}_T{\mathit{\epsilon }}_i(\mathrm{0},\mathrm{1}),$$ where ${\mathit{\epsilon }}_i(\mathrm{0},\mathrm{1})$ represents independent and normally distributed random variables with a mean of zero and a standard deviation of $\mathrm{1}$. The parameter $\mathrm{0}\le \mathit{\xi }\le \mathrm{1}$ determines the fraction of the input temperature time series which is aliased into white noise.

The covariance between the target site temperature and a precipitation-weighted temperature time series is then A8$$\mathrm{cov}\left(T_{\mathrm{0}},T_i^{\left(\mathrm{pw}\right)}\right)={\left(\mathrm{1}-\mathit{\xi }\right)}^{\mathrm{1}/\mathrm{2}}{\mathit{\sigma }}_T^{\mathrm{2}}\exp \left(-{\displaystyle \frac{r_i}{\mathit{\tau }}}\right),$$ which implies that the spatial correlation structure between $T_{\mathrm{0}}$ and the precipitation-weighted temperature follows the same exponential decay as in Eq. (), only scaled by the factor $(\mathrm{1}-\mathit{\xi }{)}^{\mathrm{1}/\mathrm{2}}$ (see dashed black line in Fig. ). The factor $\mathit{\xi }$ can be estimated from the climate model data by analysing the local correlation, i.e. at the same grid cell, between the temperature and the precipitation-weighted temperature.

We further assume that the effect of precipitation intermittency is not independent between sites but is related to the spatial coherence of the precipitation fields, for which we assume an exponential decorrelation structure with a decay length ${\mathit{\tau }}_{\mathrm{pw}}$. Based on these assumptions, the spatial covariance between sites of the white noise terms induced by the effect of precipitation intermittency has the form (see purple line in Fig. ) A9$$\mathrm{cov}\left({\mathit{\epsilon }}_i,{\mathit{\epsilon }}_j\right)=\exp \left(-{\displaystyle \frac{d_{ij}}{{\mathit{\tau }}_{\mathrm{pw}}}}\right).$$ Then, the correlation between the target site temperature and the spatial average of $N_{\mathrm{\ell }}$ precipitation-weighted temperature time series is governed by the intermittency factor $\mathit{\xi }$ and by the two spatial length scales $\mathit{\tau }$ and ${\mathit{\tau }}_{\mathrm{pw}}$, A10$$\begin{array}{rl}& \mathrm{cor}\left(T_{\mathrm{0}},{\stackrel{\mathrm{‾}}{T}}^{\left(\mathrm{pw}\right)}\right)=\\ & {\displaystyle \frac{\sqrt{\mathrm{1}-\mathit{\xi }}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}}\exp \left(-\frac{r_i}{\mathit{\tau }}\right)}{\sqrt{N_{\mathrm{\ell }}+\mathrm{2}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}-\mathrm{1}}{\sum }_{j=i+\mathrm{1}}^{N_{\mathrm{\ell }}}g\left(d_{ij};\mathit{\tau },{\mathit{\tau }}_{\mathrm{pw}},\mathit{\xi }\right)}}},\end{array}$$ with A11$$\begin{array}{rl}g\left(d_{ij};\mathit{\tau },{\mathit{\tau }}_{\mathrm{pw}},\mathit{\xi }\right):=& \left(\mathrm{1}-\mathit{\xi }\right)\exp \left(-{\displaystyle \frac{d_{ij}}{\mathit{\tau }}}\right)\\ & +\mathit{\xi }\exp \left(-{\displaystyle \frac{d_{ij}}{{\mathit{\tau }}_{\mathrm{pw}}}}\right).\end{array}$$ An example of the correlation according to Eq. () for $N_{\mathrm{\ell }}=\mathrm{2}$ and $r_{\mathrm{1}}=\mathrm{0}$ is shown as a function of $r_{\mathrm{2}}$ in Fig. .

A4 Precipitation-weighted oxygen isotope composition

For the precipitation-weighted oxygen isotope composition field, $x\equiv {\mathit{\delta }}^{\left(\mathrm{pw}\right)}$, we assume the same effect of precipitation intermittency as for the temperature field. Furthermore, an analysis of the climate model data suggests that the oxygen isotope field largely exhibits an exponential decorrelation structure in space (not shown). Hence, the correlation between the target site temperature and the spatial average of $N_{\mathrm{\ell }}$ ${\mathit{\delta }}^{\left(\mathrm{pw}\right)}$ time series is obtained in a similar manner as for $T^{\left(\mathrm{pw}\right)}$, i.e.

A12$$\begin{array}{rl}& \mathrm{cor}\left(T_{\mathrm{0}},{\stackrel{\mathrm{‾}}{\mathit{\delta }}}^{\left(\mathrm{pw}\right)}\right)=\\ & {\displaystyle \frac{\sqrt{\mathrm{1}-\mathit{\xi }}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}}\mathrm{cor}\left(T_{\mathrm{0}},{\mathit{\delta }}_i\right)}{\sqrt{N_{\mathrm{\ell }}+\mathrm{2}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{\ell }}-\mathrm{1}}{\sum }_{j=i+\mathrm{1}}^{N_{\mathrm{\ell }}}g\left(d_{ij};{\mathit{\tau }}_{\mathit{\delta }},{\mathit{\tau }}_{\mathrm{pw}},\mathit{\xi }\right)}}},\end{array}$$ where ${\mathit{\tau }}_{\mathit{\delta }}$ is the decorrelation length of the $\mathit{\delta }$ field and the only difference to Eq. () is the unknown spatial correlation structure between the temperature at the target site and the oxygen isotope field, $\mathrm{cor}(T_{\mathrm{0}},{\mathit{\delta }}_i)$. Based on our climate model results (Fig. ), we parameterize this function with a simple linear decay of the form A13$$\mathrm{cor}\left(T_{\mathrm{0}},{\mathit{\delta }}_i\right)=\left\{\begin{array}{ll}{c}_{\mathrm{0}}-\mathit{\gamma }d,& d\le d_{\mathrm{0}},\\ \mathrm{0},& d>d_{\mathrm{0}},\end{array}\right.$$ where $\mathit{\gamma }=c_{\mathrm{0}}/d_{\mathrm{0}}$, and $d_{\mathrm{0}}$ is some threshold distance above which the correlation is zero.

A5 Model parameter estimation and model results

Overall, our model is governed by three decorrelation lengths ($\mathit{\tau }$, ${\mathit{\tau }}_{\mathit{\delta }}$, ${\mathit{\tau }}_{\mathrm{pw}}$), the intermittency factor $\mathit{\xi }$, and two parameters describing the temperature–isotope correlation ($c_{\mathrm{0}}$, $d_{\mathrm{0}}$).

We estimate $\mathit{\tau }$ from the climate model data for the DML and Vostok regions (Fig. ) and find for both regions values of $\mathit{\tau }=\mathrm{1900}$ km. In the same way we estimate a value of ${\mathit{\tau }}_{\mathit{\delta }}=\mathrm{1100}$ km for both regions. The intermittency factor $\mathit{\xi }$ is derived from the local correlation between temperature and precipitation-weighted temperature (Eq. ). We find an average value for the DML region of ${\mathit{\xi }}_{\mathrm{DML}}=\mathrm{0.73}$, which is close to the average value across all of Antarctica (${\mathit{\xi }}_{\mathrm{Ant}.}=\mathrm{0.71}$), while the intermittency is stronger for the Vostok region (${\mathit{\xi }}_{\mathrm{Vostok}}=\mathrm{0.82}$). We parameterize the temperature–isotope correlation in the DML region with $c_{\mathrm{0}}=\mathrm{0.4}$ and $d_{\mathrm{0}}=\mathrm{6000}$ km and in the Vostok region with $c_{\mathrm{0}}=\mathrm{0.6}$ and $d_{\mathrm{0}}=\mathrm{2500}$ km (Fig. ). The only unconstrained parameter is the decorrelation length of the effect of precipitation intermittency, ${\mathit{\tau }}_{\mathrm{pw}}$, since it is unclear by which precipitation variable it is mainly governed (total annual amount, seasonal amount, or its distribution). An investigation with reanalysis data yielded scales between $\sim \mathrm{300}$ and $\mathrm{500}$ km for different precipitation variables , while our model data exhibit an average decorrelation length of $\sim \mathrm{600}$ km for the annual precipitation amount. Here, for the conceptual model we choose a value of $\mathrm{500}$ km.

Figure A2

Two-dimensional sampling correlation structures with temperature as predicted from our conceptual model using the model parameters from the DML (a–c) and Vostok (d–f) regions. Shown is the mean correlation of all possible single correlations for the average of two time series sampled from a pair of concentric rings around the target site for the fields of (a, d) $T_{\mathrm{2}\mathrm{m}}$, (b, e) $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and (c, f) ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$. Note that panels (a) and (d) are based on the same parameters and therefore identical.

[Figure omitted. See PDF]

We can test our assumption for the effect of intermittency based on using the estimated values of $\mathit{\tau }$ and $\mathit{\xi }$ to predict the spatial decorrelation between temperature and precipitation-weighted temperature (Eq. ). Indeed, this yields a comparably good fit to the data as an independent fit (root mean square deviation of $\sim \mathrm{0.03}$ between data and fit in both cases), supporting our assumption that intermittency can be parameterized by a partial conversion of the time series into white noise.

Similarly to analysing the climate model data, we now use our conceptual model to predict the two-dimensional ($N_{\mathrm{\ell }}=\mathrm{2}$) sampling correlation structures for the different model fields of $T_{\mathrm{2}\mathrm{m}}$, $T_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and ${\mathit{\delta }}^{\mathrm{18}}$O${}^{\left(\mathrm{pw}\right)}$ (Eqs. , , and ). Since our model space is continuous, we sample from locations placed on concentric rings around the target site. We either sample the two locations from the same ring or from two different rings using ring radii from $\mathrm{0}$ to $\mathrm{2000}$ km in increments of $\mathrm{10}$ km and calculate the average correlation for a specific ring combination. To obtain meaningful expectation values, we choose $\mathrm{36}$ locations distributed uniformly across each ring in steps of 10${}^{\circ }$, combine these locations one by one for each ring combination, and average across the correlations for each location pair. With the model parameters from the DML and Vostok regions we obtain the results displayed in Fig. , which are discussed and compared to the estimated results from the climate model data in the main text.

Code and data availability

The climate model data used in this study are freely available from the Zenodo database under 10.5281/zenodo.4001565 . Software to run the analyses and produce the figures is available as R code hosted in the public Git repository at https://github.com/EarthSystemDiagnostics/optimalcores (last access: 28 July 2021); a snapshot of this software code at the time of publication is archived on the Zenodo database under 10.5281/zenodo.5075439 .

Author contributions

TM and TL designed the research and developed the methodology. MW contributed by providing the climate model data and with his expertise on the modelling of precipitation isotopic composition. TM processed the model data, coded the analysis software, performed the analyses, and wrote the first version of the paper. All authors contributed to the interpretation of the results and to the preparation and revision of the final paper.

Competing interests

The authors declare that they have no conflict of interest.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We thank Jesper Sjolte (Lund University) for performing the ECHAM5/MPI-OM-wiso past1000 model simulation as well as Mathieu Casado, Raphaël Hébert, and Torben Kunz (AWI) for their helpful comments on this project and the paper. All plots and numerical analyses in this paper were carried out using the open-source software R: A Language and Environment for Statistical Computing. We are grateful to the two reviewers, Lenneke Jong and Dmitry Divine, and to István Hatvani and Zoltán Kern, whose comments helped to significantly improve the first version of this paper. Finally, we thank Nerilie Abram for editing the paper.

Financial support

This research has been supported by funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 716092) and by Helmholtz funding through the Polar Regions and Coasts in the Changing Earth System (PACES) programme of the Alfred Wegener Institute.​​​​​​​The article processing charges for this open-access publication were covered by the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (AWI).

Review statement

This paper was edited by Nerilie Abram and reviewed by Dmitry Divine and Lenneke Jong.