**Research article**
29 Jul 2021

**Research article** | 29 Jul 2021

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

Thomas Münch Martin Werner and Thomas Laepple

^{1},

^{2},

^{1,3}

**Thomas Münch et al.**Thomas Münch Martin Werner and Thomas Laepple

^{1},

^{2},

^{1,3}

^{1}Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, Research Unit Potsdam, Telegrafenberg A45, 14473 Potsdam, Germany^{2}Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, Bussestraße 24, 27570 Bremerhaven, Germany^{3}MARUM – Center for Marine Environmental Sciences and Faculty of Geosciences, University of Bremen, 28334 Bremen, Germany

^{1}Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, Research Unit Potsdam, Telegrafenberg A45, 14473 Potsdam, Germany^{2}Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, Bussestraße 24, 27570 Bremerhaven, Germany^{3}MARUM – Center for Marine Environmental Sciences and Faculty of Geosciences, University of Bremen, 28334 Bremen, Germany

**Correspondence**: Thomas Münch (thomas.muench@awi.de)

**Correspondence**: Thomas Münch (thomas.muench@awi.de)

Received: 24 Sep 2020 – Discussion started: 05 Oct 2020 – Revised: 21 Apr 2021 – Accepted: 29 Jun 2021 – Published: 29 Jul 2021

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 ∼ 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.

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 (Dansgaard, 1964; Lorius et al., 1969; Masson-Delmotte et al., 2008; Sjolte et al., 2011). Although the physical mechanisms that link local changes in temperature to the isotopic composition of precipitated snow are generally well understood (Dansgaard, 1964; Craig and Gordon, 1965; Jouzel and Merlivat, 1984) and can be modelled with general circulation models (Joussaume et al., 1984; Werner et al., 2011, 2016; Sjolte et al., 2011; Goursaud et al., 2018), 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 (Münch and Laepple, 2018).

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 (Dansgaard, 1964; Jouzel and Merlivat, 1984). 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 *≳*1000 km (Jones et al., 1997), which implies that ice cores
distributed on spatial scales below ∼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 (Schlosser et al., 2004; Sodemann et al., 2008; Birks and Edwards, 2009; Küttel et al., 2012). 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 (Schleiss and Smith, 2015). By this, the initial isotope signal is
weighted with the amount of precipitation, which introduces bias
(Steig et al., 1994; Laepple et al., 2011) and adds additional variability to the isotopic
record (Persson et al., 2011; Casado et al., 2020). 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
(Fisher et al., 1985; Münch et al., 2016; Laepple et al., 2016), which exhibits a small-scale
decorrelation scale of a few metres (Münch et al., 2016). 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 (Münch et al., 2016), but quantitative estimates for other Antarctic regions are still missing. A similarly high relative contribution is expected from precipitation intermittency (Laepple et al., 2018), which probably has a larger impact further inland than compared to coastal regions (Casado et al., 2020; Hatvani and Kern, 2017).

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 (Münch et al., 2016, 2017) with the signal-to-noise ratios estimated from analysing the records' temporal variability (Laepple et al., 2018) 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 (Laepple et al., 2018; Münch and Laepple, 2018). 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 (<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.1 Climate model data

We use data from the past-millennium simulation
(800–1999 CE; Sjolte et al., 2018) of the fully coupled ECHAM5/MPI-OM-wiso
atmosphere–ocean general circulation model equipped with stable isotope
diagnostics (Werner et al., 2016). 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^{∘} × 3.75^{∘}) and
with 19 vertical levels (Sjolte et al., 2018). 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 (Werner et al., 2011). These issues can be improved upon by
using a higher spatial resolution (Werner et al., 2011); 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 (Werner et al., 2016), and
the past-millennium simulation was used to reconstruct North Atlantic
atmospheric circulation in combination with ice-core isotope data
(Sjolte et al., 2018).

In this study, we use the 1200-year ECHAM5/MPI-OM-wiso time series of 2 m
surface air temperature (*T*_{2m}), precipitation (*p*), and oxygen
isotopic composition in precipitation (the relative abundance of oxygen-18 to
oxygen-16 isotopes, denoted as *δ*^{18}O) extracted from the total
number of 442 model grid cells that are available for the Antarctic continent
(Münch and Werner, 2020).

## 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 (Münch and Werner, 2020): (1) the 2 m temperature and oxygen isotopic composition
data averaged to an annual resolution without any weighting (denoted as
*T*_{2m} and *δ*^{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
*δ*^{18}O^{(pw)}).

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 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 *δ*^{18}O and
*δ*^{18}O^{(pw)} time series. This approach removes 443
anomalous annual values (<0.1 %), out of which 435 anomalies occur for the model year 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*_{2m} 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*_{ℓ}, of model grid cells and
vary *N*_{ℓ} and the distances of these locations relative to the target
site. To derive implications for actual ice-core studies, we use the
*δ*^{18}O^{(pw)} 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
*δ*^{18}O^{(pw)} with those obtained for
*T*_{2m}, ${T}_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and
*δ*^{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*_{2m} time series by sampling without
replacement *N*_{ℓ} grid cells that lie within a selection circle of
2000 km radius around the target site and then correlating the average
*δ*^{18}O^{(pw)} time series from these *N*_{ℓ} grid
cells with the target site temperature. We perform this for different *N*_{ℓ}
and determine the optimal set of cores for each value of *N*_{ℓ} 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 10^{7}, which effectively
applies to all *N*_{ℓ}≤3, or we randomly sample 10^{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*_{ℓ}>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*_{ℓ}, especially for distances farther away from
the target site due to the radially increasing number of grid cells.

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 250 km (Fig. 1) and identify all grid cells that lie
within each of these rings. The sampling of *N*_{ℓ} grid cells is then
implemented in the following two-step process: first, we determine all possible
combinations of selecting *N*_{ℓ} 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*_{ℓ} rings (see the examples in
Fig. 1b 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*_{ℓ}=2; for *N*_{ℓ}≥3, we resort to Monte Carlo sampling of the grid cell sets due to computational reasons, for which we estimated 10^{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*_{ℓ} 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*_{ℓ}=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^{∘} S, 0^{∘} E; Fig. 1a),
consisting of 26 model grid cells. This region encompasses the site of the
deep EDML ice core (EPICA community members, 2006; Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, 2016) and >50 firn and
shallow ice cores (Altnau et al., 2015). For the Vostok region, we choose an
identical latitudinal and longitudinal coverage with respect to the Vostok
station (78.47^{∘} S, 106.83^{∘} E; Fig. 1a),
covering 30 model grid cells and encompassing the sites of the deep Vostok
and Dome C ice cores, several shallower cores (Stenni et al., 2017), and the new deep drilling project (“Little Dome C”) where an ice core extending back more than 1 million years is envisaged (Passalacqua et al., 2018).

## 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.

The temperature field over Antarctica in the climate model exhibits large-scale coherent variations (Fig. 2a) 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 ∼1500 and 2500 km, while the decorrelation lengths are notably lower, with values ranging from ∼500 to 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 *δ*^{18}O^{(pw)} with the annual
*T*_{2m} time series at the same grid cell results in generally low
correlations (mean of 0.38), which across all analysed grid cells range from ∼0.1 up to ∼0.57, with ∼60 % of the correlations ≤0.4 (Fig. 2b). The correlations are overall improved when
the ${T}_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$ time series is used instead of the
*T*_{2m} time series (mean correlation of 0.53, range ∼0.1 to
0.77; Fig. 2c) but with unaffected correlation
values mostly in the coastal regions (Fig. 2d). 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 (Casado et al., 2020).

## 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*_{2m} time series at these target sites with
the spatial fields of annual temperature and of annual *δ*^{18}O,
both unweighted and weighted by the precipitation amount
(Fig. 3).

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. 3a, e). However, for
*δ*^{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 *δ*^{18}O field is not centred on the target site but
displaced by ∼1200 km towards the southeast
(Fig. 3c, 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)}$, *δ*^{18}O, and
*δ*^{18}O^{(pw)} fields
(Fig. 3f–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)}$, *δ*^{18}O, and
*δ*^{18}O^{(pw)} 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
*δ*^{18}O^{(pw)} records from grid cells within a
2000 km circle around a target site (see Sect. 2.3.2 for details) to determine what spatial array of *N*_{ℓ} ice cores optimizes the temperature correlation with the target site.

For our specific model simulation and specifying the EDML drilling site as the target site, we already know from Fig. 3a that the optimal location for a single ice core is not the local grid cell at the target site but should be a ∼1200 km southeastward site. Choosing this more distant site increases the correlation with the target temperature from an *r* value of 0.30 for the local EDML site to a value of 0.44 (Fig. 4a). 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 2000 km selection circle (Fig. 4b–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*=0.50 for *N*_{ℓ}=3, *r*=0.52 for *N*_{ℓ}=5). We obtain
comparable results when the Vostok drilling site is specified as the target
(Fig. 4d–f). The optimal single core would be at a location
∼190 km west of Vostok (*r*=0.49 compared to the local correlation of
*r*=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*=0.60 for *N*_{ℓ}=3, *r*=0.63 for *N*_{ℓ}=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 500 and 1500 km from the respective target sites, while only about 10 % lie within 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 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. 2.3.2 and the illustration in
Fig. 1), which studies the correlation patterns only as a
function of radial distance from the target site by averaging across 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. 2.3.3) 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.

When we sample only a single location (*N*_{ℓ}=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*_{2m} field itself. The average sampling correlation
structures for *T*_{2m} across the DML and Vostok regions
(Fig. 5) can be described by an exponential decay with a length scale of ∼1900 km in both cases, consistent with the estimated spatial temperature decorrelation lengths for the individual grid cells in these regions (Fig. 2a). 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 *δ*^{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. 5a). The average correlation for the innermost ring (<250 km) is ∼0.4, but it decreases only slightly within the first ∼800 km, followed by a slightly steeper decrease and nearly constant levels of *r**≲*0.2 for distances *≳*1600 km. For the Vostok region (Fig. 5b), the average sampling correlation structure
for *δ*^{18}O exhibits a nearly linear decrease from an initial
value of *r*∼0.6 to *r*∼0.1 in the final ring (>2000 km). When we
analyse the *δ*^{18}O^{(pw)} 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*_{ℓ}=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. 6). The difference in the average sampling
correlation structure between the fields of *T*_{2m} and
*δ*^{18}O^{(pw)} is even more pronounced for *N*_{ℓ}=2
than for *N*_{ℓ}=1. The maximum average correlation for *T*_{2m} is still found when both sampling locations lie inside the innermost ring, as shown for the DML region (Fig. 6a). However, for
*δ*^{18}O^{(pw)} 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 ∼1000 and 1250 km from the
target site (Fig. 6c). 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. 6b). Here, in contrast to
*T*_{2m}, 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.

Analysing the Vostok study region leads to comparable results
(Fig. 6d–f), with a similar difference in average
sampling correlation structure between *T*_{2m} 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 *δ*^{18}O^{(pw)} for distances *≲*1000 km. However, the results for the *δ*^{18}O^{(pw)} field
(Fig. 6f) 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 *δ*^{18}O field
(Fig. 5) have an influence here.

The general feature of the optimal *δ*^{18}O^{(pw)}
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 ∼82 *%* of all cases the optimal second core location is at least the second ring (>250 km), and in ∼63 % of the cases it is the second to fourth ring (250–1000 km).

We also obtain similar results when averaging *N*_{ℓ}=3 or 5 locations of
the *δ*^{18}O^{(pw)} field to reconstruct the target site
temperature (Fig. 7). 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 ∼750 and
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 (250–750 km) ring.

In summary, averaging the *δ*^{18}O^{(pw)} 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. 8a). While the local correlation for the EDML
target site stays constant around 0.31, the optimal correlation rises to
0.35 for *N*_{ℓ}=2 and to 0.43 for *N*_{ℓ}=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*_{ℓ}=2 from 0.45 to ∼0.50, but for larger *N*_{ℓ} the optimal correlation also increases more strongly and reaches 0.58 for *N*_{ℓ}=10, which is a ∼1.7-fold higher explained variance compared to *N*_{ℓ}=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*_{ℓ}=3 are actually larger than the respective local correlation (Fig. 8b).

## 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. 2b), this correlation can be increased considerably by averaging isotope records across space (Fig. 8a) following a distinct radial pattern which combines the local target site with locations between a few hundred kilometres and ∼1000 km from the target site (Figs. 6c, f, and 7). 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. 3). 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*_{2m}, ${T}_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$, and
*δ*^{18}O^{(pw)} separately in order to understand the
influence of each process (see Appendix A 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 A2). Such an exponential temperature
decorrelation is a commonly observed feature (Jones et al., 1997) and also
confirmed by our climate model data (Figs. 2a,
3a, e, and 5). 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. A2a). 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*_{ℓ}) 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 (Laepple et al., 2018; Casado et al., 2020). 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 A3). 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. A1 and A2b, 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. 6), 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 *δ*^{18}O^{(pw)} 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. 5 and Appendix A4). 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. A2c and 6c). In
addition, when we change the parameterized isotope–temperature relationship
such that it more closely resembles the Vostok region data
(Fig. 5b), the sampling correlation structure in the conceptual model (Fig. A2f) is more similar to the observed correlation structure (Fig. 6f). However, in general the conceptual model fails for
*δ*^{18}O^{(pw)} 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 *δ*^{18}O^{(pw)} 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. 2a) 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*_{ℓ}=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. 4). 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. 6) 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 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 Münch et al., 2016, 2017, 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*_{ℓ}=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 (∼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. Steen-Larsen et al., 2014; Madsen et al., 2019). 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.

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 (>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.

## 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 (Laepple et al., 2018).

In the model, we consider a temperature time series *T*_{0} at some target site
*r*_{0} and a scalar field *x* of a given climate variable. From this field, we
select *N*_{ℓ} time series *x*_{i} at the locations *r*_{i},
$i=\mathrm{1},\mathrm{\dots},{N}_{\mathrm{\ell}}$ and denote the spatial average of these time series by
$\stackrel{\mathrm{\u203e}}{x}=\frac{\mathrm{1}}{{N}_{\mathrm{\ell}}}{\sum}_{i=\mathrm{1}}^{{N}_{\mathrm{\ell}}}{x}_{i}$. The distances of
the *N*_{ℓ} 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*_{0} and $\stackrel{\mathrm{\u203e}}{x}$ follows from

and it is governed by the covariance between the temperature at the target site and the climate field at the sampling locations *r*_{i},

as well as by the covariance between the sampling locations through the variance of their spatial average,

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. A1.

## A2 Temperature

For the near-surface temperature field, *x*≡*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 *τ*; i.e. the covariance between sites is (see black line in
Fig. A1)

The correlation between the target site temperature and the spatial average of *N*_{ℓ} temperature time series is then obtained from

## A3 Precipitation-weighted temperature

To model the effect of precipitation intermittency, we follow
Laepple et al. (2018) 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 *r*_{i} arises from *T*_{i} as

where *ε*_{i}(0,1) represents independent and normally distributed random
variables with a mean of zero and a standard deviation of 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

which implies that the spatial correlation structure between *T*_{0} and the
precipitation-weighted temperature follows the same exponential decay as in
Eq. (A4), only scaled by the factor $(\mathrm{1}-\mathit{\xi}{)}^{\mathrm{1}/\mathrm{2}}$ (see dashed
black line in Fig. A1). The factor *ξ*
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 *τ*_{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. A1)

Then, the correlation between the target site temperature and the spatial
average of *N*_{ℓ} precipitation-weighted temperature time series is governed by the
intermittency factor *ξ* and by the two spatial length scales *τ* and
*τ*_{pw},

with

An example of the correlation according to Eq. (A10) for *N*_{ℓ}=2 and *r*_{1}=0 is shown as a function of *r*_{2} in
Fig. A1.

## A4 Precipitation-weighted oxygen isotope composition

For the precipitation-weighted oxygen isotope composition field, *x*≡*δ*^{(pw)}, 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*_{ℓ} *δ*^{(pw)} time series is obtained in a similar manner as for *T*^{(pw)}, i.e.

where *τ*_{δ} is the decorrelation length of the *δ* field and the only difference to Eq. (A10) is the unknown spatial
correlation structure between the temperature at the target site and the oxygen isotope field, cor(*T*_{0},*δ*_{i}). Based on our climate
model results (Fig. 5), we parameterize this
function with a simple linear decay of the form

where $\mathit{\gamma}={c}_{\mathrm{0}}/{d}_{\mathrm{0}}$, and *d*_{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 (*τ*,
*τ*_{δ}, *τ*_{pw}), the intermittency factor *ξ*, and two parameters describing the temperature–isotope correlation (*c*_{0}, *d*_{0}).

We estimate *τ* from the climate model data for the DML and Vostok regions
(Fig. 5) and find for both regions values of
*τ*=1900 km. In the same way we estimate a value of *τ*_{δ}=1100 km for both regions. The intermittency factor *ξ* is derived from the local correlation between temperature and precipitation-weighted temperature (Eq. A8). We find an average value for the DML region of *ξ*_{DML}=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 (*ξ*_{Vostok}=0.82). We parameterize the
temperature–isotope correlation in the DML region with *c*_{0}=0.4 and
*d*_{0}=6000 km and in the Vostok region with *c*_{0}=0.6 and *d*_{0}=2500 km
(Fig. 5). The only unconstrained parameter is the
decorrelation length of the effect of precipitation intermittency,
*τ*_{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 ∼300 and 500 km for different precipitation variables (Münch and Laepple, 2018), while our model data
exhibit an average decorrelation length of ∼600 km for the annual
precipitation amount. Here, for the conceptual model we choose a value of
500 km.

We can test our assumption for the effect of intermittency based on using the
estimated values of *τ* and *ξ* to predict the spatial decorrelation
between temperature and precipitation-weighted temperature
(Eq. A8). Indeed, this yields a comparably good fit to
the data as an independent fit (root mean square deviation of ∼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*_{ℓ}=2) sampling correlation structures for the different model fields of *T*_{2m}, ${T}_{\mathrm{2}\mathrm{m}}^{\left(\mathrm{pw}\right)}$,
and *δ*^{18}O^{(pw)} (Eqs. A6,
A10, and A12). 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 0 to 2000 km in increments of
10 km and calculate the average correlation for a specific ring
combination. To obtain meaningful expectation values, we choose 36 locations
distributed uniformly across each ring in steps of 10^{∘}, 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. A2, which are discussed and compared to the
estimated results from the climate model data in the main text.

The climate model data used in this study are freely available from the Zenodo database under https://doi.org/10.5281/zenodo.4001565 (Münch and Werner, 2020). 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 https://doi.org/10.5281/zenodo.5075439 (Münch, 2021).

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.

The authors declare that they have no conflict of interest.

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

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.

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).

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

Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung: Neumayer III and Kohnen Station in Antarctica operated by the Alfred Wegener Institute, Journal of large-scale research facilities, 2, A85, https://doi.org/10.17815/jlsrf-2-152, 2016. a

Altnau, S., Schlosser, E., Isaksson, E., and Divine, D.: Climatic signals from 76 shallow firn cores in Dronning Maud Land, East Antarctica, The Cryosphere, 9, 925–944, https://doi.org/10.5194/tc-9-925-2015, 2015. a

Birks, S. J. and Edwards, T. W. D.: Atmospheric circulation controls on precipitation isotope–climate relations in western Canada, Tellus B, 61, 566–576, https://doi.org/10.1111/j.1600-0889.2009.00423.x, 2009. a

Casado, M., Münch, T., and Laepple, T.: Climatic information archived in ice cores: impact of intermittency and diffusion on the recorded isotopic signal in Antarctica, Clim. Past, 16, 1581–1598, https://doi.org/10.5194/cp-16-1581-2020, 2020. a, b, c, d

Craig, H. and Gordon, L. I.: Deuterium and oxygen 18 variations in the ocean and the marine atmosphere, in: Stable Isotopes in Oceanographic Studies and Paleotemperatures, edited by: Tongiorgi, E., Proceedings Spoleto 1965, V. Lishi e F., Pisa, Italy, 9–130, 1965. a

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468, https://doi.org/10.3402/tellusa.v16i4.8993, 1964. a, b, c

EPICA community members: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198, https://doi.org/10.1038/nature05301, 2006. a

Fisher, D. A., Reeh, N., and Clausen, H. B.: Stratigraphic Noise in Time Series Derived from Ice Cores, Ann. Glaciol., 7, 76–83, https://doi.org/10.1017/S0260305500005942, 1985. a

Goursaud, S., Masson-Delmotte, V., Favier, V., Orsi, A., and Werner, M.: Water stable isotope spatio-temporal variability in Antarctica in 1960–2013: observations and simulations from the ECHAM5-wiso atmospheric general circulation model, Clim. Past, 14, 923–946, https://doi.org/10.5194/cp-14-923-2018, 2018. a

Hatvani, I. G. and Kern, Z.: Weighting alternatives for water stable isotopes: Statistical comparison between station- and firn/ice-records, Pol. Polar Res., 38, 105–124, https://doi.org/10.1515/popore-2017-0006, 2017. a

Jones, P. D., Osborn, T. J., and Briffa, K. R.: Estimating Sampling Errors in Large-Scale Temperature Averages, J. Clim., 10, 2548–2568, https://doi.org/10.1175/1520-0442(1997)010<2548:ESEILS>2.0.CO;2, 1997. a, b

Joussaume, S., Sadourny, R., and Jouzel, J.: A general circulation model of water isotope cycles in the atmosphere, Nature, 311, 24–29, https://doi.org/10.1038/311024a0, 1984. a

Jouzel, J. and Merlivat, L.: Deuterium and Oxygen 18 in Precipitation: Modeling of the Isotopic Effects During Snow Formation, J. Geophys. Res., 89, 11749–11757, https://doi.org/10.1029/JD089iD07p11749, 1984. a, b

Küttel, M., Steig, E. J., Ding, Q., Monaghan, A. J., and Battisti, D. S.: Seasonal climate information preserved in West Antarctic ice core water isotopes: relationships to temperature, large-scale circulation, and sea ice, Clim. Dyn., 39, 1841–1857, https://doi.org/10.1007/s00382-012-1460-7, 2012. a

Laepple, T., Werner, M., and Lohmann, G.: Synchronicity of Antarctic temperatures and local solar insolation on orbital timescales, Nature, 471, 91–94, https://doi.org/10.1038/nature09825, 2011. a

Laepple, T., Hörhold, M., Münch, T., Freitag, J., Wegner, A., and Kipfstuhl, S.: Layering of surface snow and firn at Kohnen Station, Antarctica: Noise or seasonal signal?, J. Geophys. Res.-Earth Surf., 121, 1849–1860, https://doi.org/10.1002/2016JF003919, 2016. a

Laepple, T., Münch, T., Casado, M., Hoerhold, M., Landais, A., and Kipfstuhl, S.: On the similarity and apparent cycles of isotopic variations in East Antarctic snow pits, The Cryosphere, 12, 169–187, https://doi.org/10.5194/tc-12-169-2018, 2018. a, b, c, d, e, f

Lorius, C., Merlivat, L., and Hagemann, R.: Variation in the Mean Deuterium Content of Precipitations in Antarctica, J. Geophys. Res., 74, 7027–7031, https://doi.org/10.1029/JC074i028p07027, 1969. a

Madsen, M. V., Steen-Larsen, H. C., Hörhold, M., Box, J., Berben, S. M. P., Capron, E., Faber, A.-K., Hubbard, A., Jensen, M. F., Jones, T. R., Kipfstuhl, S., Koldtoft, I., Pillar, H. R., Vaughn, B. H., Vladimirova, D., and Dahl-Jensen, D.: Evidence of Isotopic Fractionation During Vapor Exchange Between the Atmosphere and the Snow Surface in Greenland, J. Geophys. Res.-Atmos., 124, 2932–2945, https://doi.org/10.1029/2018JD029619, 2019. a

Masson-Delmotte, V., Hou, S., Ekaykin, A., Jouzel, J., Aristarain, A., Bernardo, R. T., Bromwich, D., Cattani, O., Delmotte, M., Falourd, S., Frezzotti, M., Gallée, H., Genoni, L., Isaksson, E., Landais, A., Helsen, M. M., Hoffmann, G., Lopez, J., Morgan, V., Motoyama, H., Noone, D., Oerter, H., Petit, J. R., Royer, A., Uemura, R., Schmidt, G. A., Schlosser, E., Simões, J. C., Steig, E. J., Stenni, B., Stievenard, M., van den Broeke, M. R., van de Wal, R. S. W., van de Berg, W. J., Vimeux, F., and White, J. W. C.: A Review of Antarctic Surface Snow Isotopic Composition: Observations, Atmospheric Circulation, and Isotopic Modeling, J. Climate, 21, 3359–3387, https://doi.org/10.1175/2007JCLI2139.1, 2008. a

Münch, T.: optimalcores: An R software project to analyse optimal ice core locations in a climate model simulation, v1.0.0, Zenodo [code], https://doi.org/10.5281/zenodo.5075439, 2021. a

Münch, T. and Laepple, T.: What climate signal is contained in decadal- to centennial-scale isotope variations from Antarctic ice cores?, Clim. Past, 14, 2053–2070, https://doi.org/10.5194/cp-14-2053-2018, 2018. a, b, c

Münch, T. and Werner, M.: Antarctic time series of temperature, precipitation, and stable isotopes in precipitation from the ECHAM5/MPI-OM-wiso past1000 climate model simulation, v0.1.0, Zenodo [data set], https://doi.org/10.5281/zenodo.4001565, 2020. a, b, c

Münch, T., Kipfstuhl, S., Freitag, J., Meyer, H., and Laepple, T.: Regional climate signal vs. local noise: a two-dimensional view of water isotopes in Antarctic firn at Kohnen Station, Dronning Maud Land, Clim. Past, 12, 1565–1581, https://doi.org/10.5194/cp-12-1565-2016, 2016. a, b, c, d, e

Münch, T., Kipfstuhl, S., Freitag, J., Meyer, H., and Laepple, T.: Constraints on post-depositional isotope modifications in East Antarctic firn from analysing temporal changes of isotope profiles, The Cryosphere, 11, 2175–2188, https://doi.org/10.5194/tc-11-2175-2017, 2017. a, b

Passalacqua, O., Cavitte, M., Gagliardini, O., Gillet-Chaulet, F., Parrenin, F., Ritz, C., and Young, D.: Brief communication: Candidate sites of 1.5 Myr old ice 37 km southwest of the Dome C summit, East Antarctica, The Cryosphere, 12, 2167–2174, https://doi.org/10.5194/tc-12-2167-2018, 2018. a

Persson, A., Langen, P. L., Ditlevsen, P., and Vinther, B. M.: The influence of precipitation weighting on interannual variability of stable water isotopes in Greenland, J. Geophys. Res., 116, D20120, https://doi.org/10.1029/2010JD015517, 2011. a

Schleiss, M. and Smith, J. A.: Two Simple Metrics for Quantifying Rainfall Intermittency: The Burstiness and Memory of Interamount Times, J. Hydrometeor., 17, 421–436, https://doi.org/10.1175/JHM-D-15-0078.1, 2015. a

Schlosser, E., Reijmer, C., Oerter, H., and Graf, W.: The influence of
precipitation origin on the *δ*^{18}O–*T* relationship at Neumayer
station, Ekströmisen, Antarctica, Ann. Glaciol., 39, 41–48,
https://doi.org/10.3189/172756404781814276, 2004. a

Sjolte, J., Hoffmann, G., Johnsen, S. J., Vinther, B. M., Masson-Delmotte, V., and Sturm, C.: Modeling the water isotopes in Greenland precipitation 1959–2001 with the meso-scale model REMO-iso, J. Geophys. Res., 116, D18105, https://doi.org/10.1029/2010JD015287, 2011. a, b

Sjolte, J., Sturm, C., Adolphi, F., Vinther, B. M., Werner, M., Lohmann, G., and Muscheler, R.: Solar and volcanic forcing of North Atlantic climate inferred from a process-based reconstruction, Clim. Past, 14, 1179–1194, https://doi.org/10.5194/cp-14-1179-2018, 2018. a, b, c

Sodemann, H., Masson-Delmotte, V., Schwierz, C., Vinther, B. M., and Wernli, H.: Interannual variability of Greenland winter precipitation sources: 2. Effects of North Atlantic Oscillation variability on stable isotopes in precipitation, J. Geophys. Res., 113, D12111, https://doi.org/10.1029/2007JD009416, 2008. a

Steen-Larsen, H. C., Masson-Delmotte, V., Hirabayashi, M., Winkler, R., Satow, K., Prié, F., Bayou, N., Brun, E., Cuffey, K. M., Dahl-Jensen, D., Dumont, M., Guillevic, M., Kipfstuhl, S., Landais, A., Popp, T., Risi, C., Steffen, K., Stenni, B., and Sveinbjörnsdottír, A. E.: What controls the isotopic composition of Greenland surface snow?, Clim. Past, 10, 377–392, https://doi.org/10.5194/cp-10-377-2014, 2014. a

Steig, E. J., Grootes, P. M., and Stuiver, M.: Seasonal Precipitation Timing and Ice Core Records, Science, 266, 1885–1886, https://doi.org/10.1126/science.266.5192.1885, 1994. a

Stenni, B., Curran, M. A. J., Abram, N. J., Orsi, A., Goursaud, S., Masson-Delmotte, V., Neukom, R., Goosse, H., Divine, D., van Ommen, T., Steig, E. J., Dixon, D. A., Thomas, E. R., Bertler, N. A. N., Isaksson, E., Ekaykin, A., Werner, M., and Frezzotti, M.: Antarctic climate variability on regional and continental scales over the last 2000 years, Clim. Past, 13, 1609–1634, https://doi.org/10.5194/cp-13-1609-2017, 2017. a

Werner, M., Langebroek, P. M., Carlsen, T., Herold, M., and Lohmann, G.: Stable water isotopes in the ECHAM5 general circulation model: Toward high-resolution isotope modeling on a global scale, J. Geophys. Res., 116, D15109, https://doi.org/10.1029/2011JD015681, 2011. a, b, c

Werner, M., Haese, B., Xu, X., Zhang, X., Butzin, M., and Lohmann, G.: Glacial–interglacial changes in H218O, HDO and deuterium excess – results from the fully coupled ECHAM5/MPI-OM Earth system model, Geosci. Model Dev., 9, 647–670, https://doi.org/10.5194/gmd-9-647-2016, 2016. a, b, c