Articles | Volume 16, issue 4
Research article
28 Jul 2020
Research article |  | 28 Jul 2020

Greenland temperature and precipitation over the last 20 000 years using data assimilation

Jessica A. Badgeley, Eric J. Steig, Gregory J. Hakim, and Tyler J. Fudge

Reconstructions of past temperature and precipitation are fundamental to modeling the Greenland Ice Sheet and assessing its sensitivity to climate. Paleoclimate information is sourced from proxy records and climate-model simulations; however, the former are spatially incomplete while the latter are sensitive to model dynamics and boundary conditions. Efforts to combine these sources of information to reconstruct spatial patterns of Greenland climate over glacial–interglacial cycles have been limited by assumptions of fixed spatial patterns and a restricted use of proxy data. We avoid these limitations by using paleoclimate data assimilation to create independent reconstructions of mean-annual temperature and precipitation for the last 20 000 years. Our method uses oxygen isotope ratios of ice and accumulation rates from long ice-core records and extends this information to all locations across Greenland using spatial relationships derived from a transient climate-model simulation. Standard evaluation metrics for this method show that our results capture climate at locations without ice-core records. Our results differ from previous work in the reconstructed spatial pattern of temperature change during abrupt climate transitions; this indicates a need for additional proxy data and additional transient climate-model simulations. We investigate the relationship between precipitation and temperature, finding that it is frequency dependent and spatially variable, suggesting that thermodynamic scaling methods commonly used in ice-sheet modeling are overly simplistic. Our results demonstrate that paleoclimate data assimilation is a useful tool for reconstructing the spatial and temporal patterns of past climate on timescales relevant to ice sheets.

Please read the corrigendum first before continuing.

1 Introduction

Predicting the future behavior of the Greenland Ice Sheet requires understanding its sensitivity to changes in temperature and precipitation (Bindschadler et al.2013). One important constraint on this sensitivity is the response of the paleo ice sheet to changing climate in the past. On glacial–interglacial timescales, temperature, not precipitation, appears to be the dominant control on the size of the Greenland Ice Sheet (Alley et al.2010), as evidenced by the fact that the ice sheet is largest during cold and arid glacial periods and smallest during warm and wet interglacials. On these timescales, precipitation over the Greenland Ice Sheet scales positively with temperature (Robin1977), as anticipated by the Clausius–Clapeyron relation between temperature and saturation vapor pressure. Ice-core records, however, show that this thermodynamic relation is a poor approximation on annual to multi-millennial timescales (Kapsner et al.1995; Fudge et al.2016). For example, the GISP2 ice core from central Greenland shows that cooling coincided with increased snowfall between the early Holocene and present (Cuffey and Clow1997). Despite such evidence, paleo ice-sheet models typically assume precipitation fields that are parameterized in time using a thermodynamic relationship that is constant for all locations and timescales (e.g., Huybrechts2002; Greve et al.2011).

Ice-core records provide the best empirical estimates of climate history over the Greenland Ice Sheet. For temperature, important proxies include oxygen and hydrogen isotopes of ice (e.g., Jouzel et al.1997), nitrogen isotope ratios of gas trapped in ice (e.g., Severinghaus et al.1998; Severinghaus and Brook1999), and borehole thermometry (e.g., Cuffey et al.1995; Dahl-Jensen et al.1998). For precipitation, the thickness of annual layers of accumulated ice, corrected for thinning, is used (e.g., Dahl-Jensen et al.1993). Ice-core records, however, are limited in their spatial coverage. In contrast, climate-model simulations are spatially complete estimates of past climate, but they are subject to uncertainty due to model dynamics, boundary conditions, and spatial resolution.

Efforts to combine information from proxy data and climate models have long been a part of ice-sheet modeling. The most common approach is to scale the modern spatial pattern of temperature and precipitation using data from a single ice core (e.g, Huybrechts et al.1991; Huybrechts2002; Greve1997; Greve et al.2011; Nielsen et al.2018). This assumes a fixed spatial pattern through time, which is unlikely to be valid. Recently, Buizert et al. (2018) used the average of the three best-understood Greenland ice-core records to adjust the results of a transient climate-model simulation (the transient climate evolution experiment, TraCE-21ka; Liu et al.2009; He et al.2013). This approach allows for possible changes in spatial relationships but focuses on ice cores in central and northern Greenland and provides no information on precipitation. Other attempts to incorporate more proxy data have been limited to short time periods (e.g., Simpson et al.2009; Lecavalier et al.2014).

In this study we apply paleoclimate data assimilation to obtain a new, spatially complete reconstruction of Greenland temperature and precipitation. We focus on the last 20 000 years, which include the end of the last glacial period, the glacial to interglacial transition, and the Holocene Thermal Maximum (HTM), when temperatures at the Greenland Ice Sheet summit reached 1–2 C warmer than the present (Cuffey and Clow1997; Dahl-Jensen et al.1998). The climate history and the corresponding changes in the size of the ice sheet are well-documented over this time period (e.g., Kaufman et al.2004; Young and Briner2015).

Paleoclimate data assimilation combines spatial information from a climate-model simulation and temporal data from proxy records to produce a climate “reanalysis”, where the term is taken from the modern climate reanalysis methods on which the data assimilation method is based (e.g., Kalnay et al.1996). We adapt the paleoclimate data assimilation framework developed by Hakim et al. (2016), who reconstructed annual 2 m air temperature and 500 hPa geopotential height over the last millennium using a global network of temperature- and precipitation-sensitive proxy records. Here, we use oxygen isotopes of ice and layer thickness from ice cores to reconstruct temperature and precipitation, respectively. We choose these proxies for their high temporal resolution, direct relationships to climate over the ice sheet, and availability from multiple ice cores. For the climate-model simulation, we use TraCE-21ka (Liu et al.2009; He et al.2013), which was run with the Community Climate System Model version 3 (CCSM3) to simulate the last 22 000 years. We compare the resulting reanalysis to previously published climate reconstructions (Sects. 3.1, 4.2, and 4.3) and assess the precipitation–temperature relationship (Sect. 4.1). We evaluate the reanalysis with independent proxy records and sensitivity tests (Sect. 3.2 and 3.3).

2 Methods

Our paleoclimate reconstruction method assimilates oxygen isotope ratios and accumulation from ice cores with a transient climate-model simulation to reconstruct the last 20 000 years of Greenland temperature and precipitation. In the following subsections we discuss the ice-core data, the climate-model simulation, and the details of our paleoclimate data assimilation approach.

2.1 Ice-core data

We use proxy records from eight ice cores from the Greenland Ice Sheet and nearby ice caps (Fig. 1, Table 1). As a proxy for temperature, we use previously published measurements of oxygen isotope ratios from the ice, which we discuss using the conventional δ18O nomenclature (Dansgaard1964). We note that while other proxies (such as borehole thermometry or δ15N of N2) have been used to produce temperature estimates (e.g., Cuffey et al.1995; Dahl-Jensen et al.1998; Severinghaus et al.1998; Severinghaus and Brook1999), they are not available for many of the core locations; we instead rely on such data to obtain independent estimates of error in the δ18O–temperature relationship (see Sect. 2.3.3) and as comparisons to our resulting reanalysis (see Sect. 3). The accumulation history has been estimated for five of these cores from layer thickness corrected for vertical ice thinning due to dynamical strain in the ice sheet. We rely on previously published accumulation histories for the GISP2 and NEEM cores (Cuffey and Clow1997; Rasmussen et al.2013), and we estimate accumulation for the Dye3, GRIP, and NGRIP cores using available layer-thickness data and simple ice-thinning calculations (see below and Sect. S2 in the Supplement). We do not use accumulation records from the Agassiz, Camp Century, or Renland cores because the ice-thinning history at these sites is not adequately constrained. Most of the ice-core data are available at 50-year or higher resolution and have been synchronized to a common depth–age scale (the Greenland ice-core chronology 2005, GICC05; Andersen et al.2006; Rasmussen et al.2006; Svensson et al.2006; Vinther et al.2006, 2008). All of these ice-core records extend from the beginning of the Holocene to the present. Five δ18O and four accumulation records also include the last glacial period. To evaluate the impact of the differing lengths of these records, we produce a sensitivity reanalysis for which we assimilate just the fixed proxy network (i.e., only those data that span the full reanalysis time period, the last 20 000 years).

Figure 1Locations of the ice-core sites referenced in this study. We use oxygen isotope (δ18O) records from all eight sites and accumulation records from the five circled sites.

Table 1Metadata for the water isotope (δ18O) and accumulation (accum) records referenced in this study. “NBI” refers to the Niels Bohr Institute data access site (, last access: 24 June 2019) and “Pangaea” refers to the Pangaea data access site (, last access: 20 August 2019). Latitude and longitude are in units of decimal degrees (dd) and dates are in thousands of years before 1950 CE (ka).

Download Print Version | Download XLSX

To extract the accumulation signal from layer thickness, the layers must be destrained using assumptions about the history of ice flow. For the Dye3, GRIP, and NGRIP cores, we use a one-dimensional ice-flow model (Dansgaard and Johnsen1969) to calculate the cumulative vertical strain the layers have experienced at each core site. The Dansgaard–Johnsen model requires specifying the vertical velocity at the surface and a kink height which determines the shape of the vertical velocity profile. The velocity profile below the kink height approximates the influence of greater deformation rates in deeper ice due to increased deviatoric stress and warmer ice temperature. For sites at the pressure-melting point at the bed, such as NGRIP, we also implement the basal melt rate (e.g., Dahl-Jensen et al.2003). Previous work on many of the Greenland ice cores has estimated cumulative vertical thinning assuming that the accumulation history scales with δ18O and then found optimal parameter values by comparing the modeled and measured depth–age relationships (Dahl-Jensen et al.1993, 2003; Rasmussen et al.2014). Here, we wish to maintain independent determinations of the δ18O and accumulation proxy records. To do this, we select reasonable ice-flow parameters independently, based on the glaciological setting of each site; specifically, we use kink-height values of 0.1–0.2 for flank flow and 0.4 for ice flow near ice divides where the deviatoric stress is low (Raymond1983; Conway et al.1999). Where available, we use published values or kink-height values that result in a good match to published accumulation records (Dahl-Jensen et al.2003; Gkinis et al.2014). Based on this range of plausible ice-flow parameters, we develop three scenarios for each site: “low”, “moderate”, and “high”, where the labels reflect the relative magnitude of accumulation in the glacial and early Holocene. We assimilate the intermediate-value (moderate) accumulation records to produce our main precipitation reanalysis, while we assimilate the high- and low-accumulation records into high- and low-sensitivity scenarios, respectively, to provide a conservative estimate of uncertainty. Descriptions of the rationale for the parameter choices at each site are given in Sect. S2. Our method to estimate accumulation should be most accurate for the interior ice cores (i.e., GISP2, GRIP, NEEM, NGRIP); these sites are thicker and have lower accumulation rates such that layers of the same age have experienced less cumulative strain than for the more coastal cores (i.e., Agassiz, Camp Century, and Renland). We do not attempt to reconstruct accumulation from these coastal cores because the layers cannot be destrained with sufficient accuracy. Dye3 suffers from the same challenges; however, it is the only ice core with long-term climate data south of 70 N (Fig. 1). Thus, we include the Holocene Dye3 accumulation rates despite the greater uncertainty relative to the interior cores.

Because records from the Dye3 ice core are our only source of information in southern Greenland, we take the following steps to increase the data available for assimilation. The Dye3 record has not been previously assigned a depth–age scale beyond 11.7 ka (throughout this paper, “ka” refers to thousands of years before 1950 CE). We extend the depth–age scale to 20 ka to take advantage of the glacial portion of the δ18O record. To do this, we match the δ18O record from Dye3 to the δ18O record from NGRIP using the cross-correlation maximization procedure from Huybers and Wunsch (2004) (Sect. S3). We interpolate the glacial δ18O record from Dye3 (which, as measured, has an average resolution of only 85 years) to the same 50-year resolution used for our other ice-core records. Extension of the Dye3 depth–age scale also provides a layer-thickness record from 20 ka to the present; however, we do not use accumulation data from Dye3 for the period 20–11.7 ka because the low resolution impedes our ability to estimate accumulation variations from layer thickness. Using this depth–age scale extension for Dye3 may introduce error that is difficult to quantify; however, we find that including Dye3 has an important impact on the resulting reanalysis of southern Greenland climate (Sect. S4).

Where possible, we account for non-local effects on the ice-core records. The global-mean δ18O of seawater fluctuates with global ice volume, while on the regional scale, horizontal advection brings ice from other elevations and latitudes. We correct for changes in the oxygen isotope composition of seawater following the methods of Stenni et al. (2010), using the benthic foraminifera dataset from Bintanja et al. (2005). For ice cores in regions of high horizontal advection, we correct for elevation and latitude differences between the site of snow deposition and the ice-core site. Following the methods from Dahl-Jensen et al. (2013), we apply corrections for advection-caused elevation changes in the δ18O records from Camp Century, Dye3, and NEEM and for advection-caused latitude changes in the δ18O record from NEEM (Vinther et al.2009; Dahl-Jensen et al.2013). We do not correct the accumulation records for advection from upstream because the elevation–accumulation relationship is complicated by the prevailing wind direction (Roe and Lindzen2001) and the thinning function uncertainties are likely to be larger than the effects of ice advection. We also do not correct δ18O or accumulation for changing elevation at the ice-core site; our goal is to reconstruct conditions at the surface rather than at a constant reference elevation. We take the anomaly of each corrected δ18O record and the ratio of each accumulation record relative to the mean of all data in the record that falls within the time period 1850–2000 CE. We then average each record to 50-year resolution, the lowest resolution in these records (with the exception of δ18O from the glacial period in the Dye3 core). It is these corrected, averaged records that we use in the data assimilation (Figs. 2 and 3).

Figure 2δ18O records assimilated into the temperature reconstruction. Records are shown as anomalies relative to the mean during 1850–2000 CE and are ordered top to bottom from northernmost to southernmost. Ice-core site names are given above each record.


Figure 3Accumulation records assimilated to reconstruct precipitation for the main reanalysis and two sensitivity scenarios. Records are shown as fractions relative to the mean during 1850–2000 CE and are ordered top to bottom from northernmost to southernmost. Black lines are the moderate records which are included in the main precipitation reanalysis, red lines are the low records which are included in the low-sensitivity scenario, and blue lines are the high records which are included in the high-sensitivity scenario. Note that we use the same GISP2 accumulation record for the main, high, and low scenarios. Ice-core site names are given above each set of records.


2.2 Climate-model simulation

We use TraCE-21ka, a simulation of the last 22 000 years of climate (22 to −0.04 ka), which was run using the fully coupled CCSM3 at T31 resolution (approximately 3.75 horizontally) with transient ice-sheet, orbital, greenhouse gas, and meltwater flux forcings (Liu et al.2009, 2012; He et al.2013). For paleoclimate data assimilation, it is important that the climate simulation capture a range of possible climate states over the time period of interest. By design, TraCE-21ka captures the major glacial-to-Holocene temperature changes, as well as some of the short-term, rapid climate changes, such as the Bølling–Allerød transition (Liu et al.2009). Many higher-resolution climate simulations are transient only over the last millennium (e.g., PAGES 2k-PMIP3 group2015) or provide a snapshot of a certain time, such as the Last Glacial Maximum or the mid Holocene (e.g., Harrison et al.2014). Individually, these millennial-length simulations have too little variability to capture the range of climate states across the glacial–interglacial transition. If combined, the biases in each simulation would need to be individually addressed, which is beyond the scope of this study.

From TraCE-21ka, we use 2 m air temperature for temperature (T) and the sum of large-scale stable precipitation and convective precipitation for precipitation (P). To correct for model bias in TraCE-21ka, we assume that the bias is stationary in time and apply the delta-change method (Teutschbein and Seibert2012) by taking the anomaly of temperature and the fraction of precipitation relative to the mean of our reference period (1850–2000 CE). An assumption of a stationary model bias is required because, with a small number of proxy records, we cannot afford to subsample them for the purposes of bias correction, data assimilation, and evaluation. After the bias correction, we average the TraCE-21ka variables (which originally have monthly resolution) to 50-year resolution, as we did for the ice-core records. In this process, we average 50 consecutive years (600 months) such that no year (or month) is used in more than one 50-year average. This averaging results in 440 time steps spaced 50 years apart.

TraCE-21ka includes changes in orbital forcing, which contribute to changes in the seasonality of temperature and precipitation. The strength of these seasonal cycles influences the mean-annual relationship between δ18O and temperature (Steig et al.1994; Werner et al.2000; Krinner and Werner2003). TraCE-21ka consistently shows stronger temperature and precipitation seasonality in the glacial (20 to 15 ka) than in the Holocene (5 to 0 ka) at each of the eight ice-core sites considered in this study (Fig. S3). Relative to the annual mean, the glacial had warmer and wetter summers and colder and drier winters. The Holocene also shows such a seasonal cycle; however, there is a smaller difference between the summers and winters. Any seasonal signal with wetter summers than winters will bias the δ18O towards summer values. According to TraCE-21ka, this effect is amplified in the glacial. As we discuss in Sect. 2.3.1, the particularly strong summer bias in the glacial affects the mean-annual δ18O–temperature relationship in ways that are consistent with findings from borehole thermometry at the GISP2, GRIP, and Dye3 ice-core sites (Cuffey and Clow1997; Jouzel et al.1997).

In addition to a change in the strength of the seasonal cycle, TraCE-21ka shows a temporal shift, with summer temperature and precipitation peaking earlier in the glacial (around June and July) than in the Holocene (from July to September) (Fig. S3). In the glacial, both variables peak around June and July, with only two exceptions: precipitation peaks in August at the Renland and Dye3 ice-core sites. In contrast, Holocene temperature peaks slightly later, in July, while precipitation peaks even later, in August and occasionally September. Both variables and both time periods show winter minimums in February, again with the two exceptions of Renland and Dye3, which show later precipitation minimums. In this study, the timing of the seasonal peaks is relevant because it affects the precipitation-weighted temperature, defined in Eq. (7).

TraCE-21ka also includes prescribed transient ice sheets as a boundary condition. The transience is important for capturing the influence of elevation change on the ice-core records; however, the low horizontal resolution of TraCE-21ka leads to difficulties in capturing elevation changes at the edges of the ice sheet, in southern Greenland, and at coastal ice caps. In addition, the ice sheets in TraCE-21ka are independent of climate, updated only every 500 years during the simulation, and taken from ICE-5G (Peltier2004), a now outdated ice-sheet reconstruction (Roy and Peltier2018).

2.3 Paleoclimate data assimilation

To combine the ice-core and climate-model data, we use an offline data assimilation method similar to that described in Hakim et al. (2016). If no covariance localization is used, as in this study, this method can be summed up as a linear combination of randomly selected model states that are weighted according to new information provided by the proxy records. “Offline” refers to the absence of a forecast model that evolves the climate state between assimilation time steps. The offline method is appropriate when model predictive skill is small given the assimilation time step (Hakim et al.2016, and references therein). Model predictive skill is generally poor on decadal to longer timescales (Latif and Keenlyside2011, and references therein) except possibly during times of strong forcing, such as the Bølling–Allerød (14.7–12.7 ka) and the Younger Dryas (12.7–11.7 ka) (Hawkins and Sutton2009). Because each of our time steps is an average over 50 years, as dictated by the resolution to which we average the proxy records, the offline method is appropriate except possibly during these large-forcing events. For these events, an online method may be appropriate (assuming that the models correctly capture both the forcing and the response); however, online ensemble data assimilation over glacial–interglacial cycles using a fully coupled earth system model is impractical due to the computational cost and is beyond the scope of this study.

Our paleoclimate data assimilation framework uses ensembles for the initial (prior) and final (posterior) estimates of the climate state, providing a probabilistic framework for interpreting and evaluating the results. To compute the posterior ensemble, we apply the Kalman update equation (Whitaker and Hamill2002), which spreads new information gained from proxy records to all locations and variables in the prior ensemble:

(1) x a = x b + K y - H ( x b ) .

Bold lowercase letters are vectors, bold capital letters are matrices, and script capital letters are mapping functions. The posterior ensemble is xa, xb is the prior ensemble, y is the proxy data, and is the function that maps from the prior variables to the proxy variables. K is the Kalman gain matrix:

(2) K = BH T HBH T + R - 1 ,

where T indicates a matrix transpose, B is the covariance matrix computed from the prior ensemble, H is the linearization of about the mean value of the prior, and R is a diagonal matrix containing the error variance for each proxy record, the use of which requires an assumption that error covariances between proxy records are negligible.

To compute the new information gained from the proxy records, the prior ensemble must first be mapped into proxy space to get prior estimates of the proxy (ℋ(xb)). This mapping () is the proxy system model (PSM). Our PSM for the δ18O–temperature relationship is a linear function and for accumulation is a direct comparison between ice-core-derived accumulation and precipitation from the prior. For these PSMs, both temperature and precipitation are interpolated from the climate-model grid to the geographic location of the ice core. These PSMs are detailed in Sect. 2.3.1 and 2.3.2. Comparing the prior estimates of the proxy (ℋ(xb)) to the proxy data (y) yields the innovation (y−ℋ(xb)), the new information gained from the proxy records.

The Kalman gain (K, Eq. 2) weights the innovation by the relative magnitude of the ensemble covariance of the prior estimates of the proxy (HBHT) and the error covariance of the proxy records (R). The Kalman gain spreads the weighted innovation to all locations and variables in the prior ensemble, using the covariance structure (BHT) from the prior ensemble.

The prior ensemble is an initial estimate of possible climate states, which we form using 100 randomly chosen 50-year averages from the TraCE-21ka simulation. States from both the glacial and the Holocene make up a prior ensemble. The same prior is used for all time steps in the reconstruction, leading to a prior that is constant in time. For a longer discussion on the reasoning behind our choice to use a constant prior, please refer to the Supplement, Sect. S1. Proxy records are assimilated into the prior using Eq. (1), which produces the posterior ensemble, a new estimate of possible climate states. We assimilate δ18O to reconstruct temperature and separately assimilate accumulation to reconstruct precipitation. This approach maintains independence between temperature and precipitation, which avoids imposing linearity and stationarity on the relationship between these two variables. As Cuffey and Clow (1997) show, not only is this relationship non-linear on long timescales, but it is also not well-approximated by simple thermodynamic expectations. Separating these variables ensures that the relationship between temperature and precipitation is consistent with the empirical relationship between δ18O and accumulation from ice cores rather than being derived exclusively from the climate model.

We repeat the data assimilation process over multiple iterations, with each iteration using 1 of 10 different 100-member prior ensembles and excluding one proxy record. Each of the 10 prior ensembles is made up of a different random selection of 50-year averages from TraCE-21ka. Thus, each prior ensemble has a different variance and spatial covariance structure. Each proxy record is excluded from a total of 10 iterations, where each of these iterations uses a different one of the 10 prior options. Every iteration is uniquely identifiable by which prior ensemble is used and which proxy record is excluded. For a reanalysis, the total number of iterations is thus 10 times the number of proxy records such that for temperature we have 80 iterations and for precipitation we have 50 iterations. A reanalysis is a compilation of the 100-member posterior ensembles from these iterations, resulting in a temperature reanalysis having 8000 ensemble members and a precipitation reanalysis having 5000 ensemble members.

The proxy record that is excluded from an iteration is independent of that iteration's posterior ensemble such that we can evaluate the posterior against this record. With our PSMs, we convert the posterior ensemble into predictions of the independent record using the mapping and compare these predictions to the record along the time axis. We use four skill metrics to evaluate different aspects of the predictions. The correlation coefficient (corr; Eq. 3) measures the relative timing of signals in the predictions and the proxy record:

(3) corr = 1 n - 1 i = 1 n ( y i - y ) ( v i - v ) σ y σ v ,

where v is the ensemble mean of the predicted values, y is the proxy record value, an overbar indicates a time mean, n is the number of time steps, and σ is the standard deviation of the variable in the subscript. The coefficient of efficiency (CE; Eq. 4) (Nash and Sutcliffe1970) is affected by signal timing as well as signal amplitude and mean bias:

(4) CE = 1 - Σ i = 1 n ( v i - y i ) 2 Σ i = 1 n ( v i - v ) 2 .

The root mean squared error (RMSE; Eq. 5) gives an intuitive sense of the magnitude of the differences between the predictions and proxy record:

(5) RMSE = 1 n i = 1 n ( y i - v i ) 2 1 / 2 .

The ensemble calibration ratio (ECR; Eq. 6) indicates whether the ensemble has enough spread (uncertainty) given the error in the ensemble mean (e.g., Houtekamer et al.2005):

(6) ECR = 1 n i = 1 n ( y i - v i ) 2 var ( v i ) + r ,

where v is a vector of the ensemble members of the predicted values, r is the error variance for the proxy record (y), and var indicates the variance. Accordingly, if the ensemble variance is appropriate for the amount of error, then the ensemble calibration ratio is near unity.

We compute all four skill metrics for both the posterior and prior ensembles, which shows whether assimilation of proxy records results in an improved estimate of the climate state over our initial estimate. We define improvement as a correlation coefficient closer to 1, CE closer to 1, RMSE closer to 0, and ensemble calibration ratio closer to 1. We anticipate that our reanalysis will show improvement over the prior because the prior is constant in time and contains no information about temporal climate variations; however, improvement is not guaranteed, especially if proxy records contain highly localized climate signals or if the prior covariance structure is unable to appropriately spread information from the proxy sites to other locations. For further comparison, we additionally compute the correlation coefficient, CE, and RMSE between the TraCE-21ka simulation and the proxy records.

We evaluate results over three time periods: (1) the full overlap between the reanalysis time period and the proxy record, (2) a time representative of the glacial (20–15 ka), and (3) a time representative of the Holocene (8–3 ka). Some proxy records overlap the full reanalysis period (20–0 ka), while others overlap just the Holocene (11.7–0 ka). The skill metrics computed for these two groups should be considered separately.

To remove mean bias from temperature, we subtract out the reference-period mean. For precipitation, we divide by the reference-period mean. It is these bias-corrected results that are referred to unless noted otherwise.

2.3.1 Proxy system model: δ18O

The isotopic composition of precipitation, as recorded in ice cores, is highly correlated with temperature at the time and location of deposition but is also sensitive to conditions at the moisture source region (i.e., sea surface temperature and relative humidity at the ocean surface; e.g., Johnsen et al.1989). Moisture source conditions primarily affect the deuterium excess, which we do not use here, and the influence on δ18O is comparatively weak (e.g., Armengaud et al.1998). For our δ18O PSM, we use a linear relationship with temperature at the ice-core drill site (Tsite) that has a slope of 0.67±0.02 ‰ C−1, which was calibrated using modern, spatial data (Johnsen et al.1989). It is well known that this modern, spatially derived slope does not necessarily apply to temporal δ18OTsite relationships, which have effective slopes that are time, frequency, and location dependent. Temporal changes in precipitation seasonality, inversion-layer thickness, and source region conditions introduce nonlinearity into the effective δ18OTsite relationship (e.g., Jouzel et al.1997; Pausata and Löfverström2015). Diffusion in the firn column also affects this relationship, but it is negligible for annual and longer timescales at the locations of the ice cores we use (Cuffey and Steig1998). Borehole thermometry at the GISP2 and GRIP sites shows that for the low-frequency changes associated with the last glacial–interglacial transition, the temporal slope is less than half the modern, spatial slope (Cuffey and Clow1997; Jouzel et al.1997). Numerous studies have suggested that precipitation seasonality is the largest source of nonlinearity in the δ18OTsite relationship (e.g., Steig et al.1994; Pausata and Löfverström2015); changes in precipitation seasonality are thought to be the primary reason that the effective δ18OTsite relationship for the glacial–interglacial transition has such a low slope (Werner et al.2000; Gierz et al.2017; Cauquoin et al.2019).

We rely on TraCE-21ka to estimate the site-specific effects of precipitation seasonality on the δ18OTsite relationship. Site-specific effects can also be estimated using independent temperature reconstructions, e.g., from borehole thermometry or δ15N of N2 measurements; however, such independent reconstructions for the last 20 000 years exist only for a few of the ice cores, GISP2, GRIP, and NGRIP (Buizert et al.2018; Dahl-Jensen et al.1998; Gkinis et al.2014), limiting the utility of such records to capture the spatial variability of the δ18OTsite relationship across Greenland.

To incorporate estimates of the site-specific effects of precipitation seasonality from TraCE-21ka, we adjust the linear δ18O PSM by replacing Tsite with Tsite, the precipitation-weighted temperature at the model grid cell closest to the ice-core drill site. We compute T across Greenland using T and P from TraCE-21ka (see Fig. S4 for a visual comparison of T and T at the GISP2 ice-core site):

(7) T = i = 1 n = 12 T mon P mon P ann .

With Tsite in our PSM, we find that the δ18OTsite slope is spatially variable (e.g., Fig. S5), ranging between 0.42 ‰ C−1 and 0.66 ‰ C−1 at the ice-core sites (Table S1) and tending to be less than the modern spatial relationship of 0.67 ‰ C−1 at most locations around Greenland. Due to the data assimilation method outlined above, these slopes vary both in space and across iterations, the latter being due to the varying prior ensembles. These slopes do not vary in time in the prior, but they do in the posterior (note that δ18OTsite slopes mentioned throughout this paper refer to the prior ensemble). By using 10 different prior ensembles, we capture the uncertainty in the δ18O–temperature relationship due to model variations in the precipitation seasonality. These TraCE-21ka-derived estimates lie within the range of slopes estimated for sites around Greenland for a variety of time periods (Table S1). Differences seen in Table S1 reflect both the different methods used and the time period considered. Some estimates, such as Guillevic et al. (2013) and Buizert et al. (2014), are for abrupt transitions, such as Dansgaard–Oeschger events, while others find mean slopes over longer periods of time, such as Kindler et al. (2014) and this investigation.

To evaluate the sensitivity of our results to the choice of PSM, we produce four other reconstructions (S1–S4). The S1 scenario uses the PSM, δ18O=0.67Tsite, which is the modern (high-frequency) relationship and does not account for precipitation seasonality. The S2 scenario uses δ18O=0.5Tsite, the mean of the high-frequency and low-frequency temporal slopes. The S3 scenario uses δ18O=0.335Tsite, which is similar to published estimates of the glacial–interglacial temporal slope (half the high-frequency slope) (Cuffey and Clow1997; Jouzel et al.1997). By lowering the slope in the S2 and S3 scenarios, we implicitly account for precipitation seasonality; however, in these scenarios and S1, the δ18OTsite relationship is spatially uniform, whereas it is spatially variable in the main reanalysis because we include the spatial pattern of precipitation seasonality. The S4 scenario uses the same PSM as in the main reanalysis, δ18O=0.67Tsite, but we adjust the strength of the precipitation seasonality in TraCE-21ka such that the average δ18OTsite slope around Greenland is approximately 0.335 ‰ C−1. The S4 scenario thus has the same spatial variability in the δ18OTsite relationship as in the main reanalysis but a greater influence of precipitation seasonality.

These sensitivity tests are equivalent to testing different assumptions about the δ18OTsite relationship. The availability of a 20 000-year-long isotope-enabled climate simulation would allow us to determine this relationship from model physics, which incorporate a variety of processes that can affect water isotopes, including precipitation seasonality.

2.3.2 Proxy system model: accumulation

Accumulation is closely related to total precipitation at our ice-core sites, which have limited surface melting and evaporation. Simulations from the regional climate model HIRHAM5 show that for modern climate at the GISP2, GRIP, NEEM, and NGRIP ice-core sites, surface mass balance, snowfall, and precipitation are all within 1.6 cm water equivalent (w.e.) when averaged over the years 1989 to 2012 CE (Langen et al.2015, 2017). For this reason, and because TraCE-21ka lacks process-based ablation variables, our PSM is a direct comparison between ice-core accumulation and simulated precipitation at the model grid cell closest to the ice-core site.

This direct-comparison PSM may be an incomplete model of the accumulation–precipitation relationship at the Dye3 ice-core site; regional climate simulations show that modern surface mass balance is lower than precipitation due to melt rates that average 84 cm w.e. yr−1. Significant melt rates would cause the spatial covariance structure of accumulation across these sites to differ from that of precipitation; however, both models and observations lack the necessary variables or duration to show the extent to which this difference exists for 50-year timescales through the last glacial–interglacial cycle. We assimilate relative rather than absolute changes to account for the mean bias between precipitation and accumulation.

2.3.3 Proxy error variance estimation

In the Kalman filter (Eq. 1), the diagonal elements of R contain the error variance of each proxy record, which includes how we model the proxy (the PSM). We compute representative error variances for δ18O and accumulation and apply them to all records and time slices. We do not include error associated with corrections applied to the ice-core data (Sect. 2.1) or associated with the accumulation PSM (Sect. 2.3.2) because we cannot characterize the statistical properties of these errors.

A universal, but typically small, error source is from the measurement of proxies. For δ18O, measurement error is equivalent to laboratory precision. We compute a representative measurement error from the GISP2 ice core, for which a single measurement of δ18O has a laboratory precision (variance) of 0.024 ‰2 (Stuiver and Grootes2000). Assuming independent error and annual measurements, the 50-year average error variance reduces to 0.0034 ‰2, which is insignificant compared to other sources of error. For accumulation, the measurement error is from the measurement of layer thickness, which is related to the error in annual-layer counts per unit depth. Again, we assume GISP2 is a representative core and estimate the layer-thickness error from Table 3 in Alley et al. (1997), which provides repeat annual-layer counts for several depth intervals. From this table, we find the standard deviation of counted years in each depth interval, divide by the average number of years, average across all depth intervals, and square the result. This computation results in a layer-thickness error variance of 0.0015, a unitless number due to our use of fractional accumulation records.

Another source of error is the extent to which a model grid cell may misrepresent a point proxy measurement. In the innovation, there is an implicit assumption that the proxy (y) and the prior estimate of the proxy (ℋ(xb)) are representative of the same processes. However, an ice core is about 100 cm2, an area that is affected by processes at all scales, from regional change to local, sub-meter-scale topography, while a model grid cell in TraCE-21ka can cover tens of thousands of square kilometers, an area that is affected by only the largest scales, from global to regional. Thus, there is an inherent inability of the prior to represent local processes at the ice-core site, which we refer to as the spatial representation error. To estimate this error, we compute the variance of the local noise (e.g., Reeh and Fisher1983) using the GISP2 and GRIP ice-core records, which are located about 30 km apart within the same model grid cell. For δ18O, our estimate is 0.21 ‰2, which is about half the value determined by Fisher et al. (1985) at several locations around Greenland. Our estimate is relatively conservative, considering that we are using 50-year averages rather than annual averages as in Fisher et al. (1985). For accumulation, we estimate a spatial representation error variance of 0.0023 using the same method as for δ18O.

A third source of error is the extent to which the δ18O PSM may be an inaccurate model of the δ18OTsite relationship. The less accurate the PSM, the less weight should be given to the innovation. We estimate PSM error variance by calculating the mean squared error (MSE) between a δ18O record and an independent temperature record mapped to δ18O using the PSM (δ18O=0.67Tsite). We use independent datasets taken from the GISP2 ice core: the δ18O record and three δ15N-derived temperature estimates for the Holocene, including a mean estimate and the 2 standard deviation uncertainty bounds (Kobashi et al.2017). From these datasets, we estimate three PSM error variances that range from 0.56 ‰2 to 1.1 ‰2, from which we choose the largest.

For the assimilation of each δ18O record, we use an estimated total error variance of 1.3 ‰2, which is a sum of the measurement, spatial representation, and PSM error variances. For the assimilation of each accumulation record, we use an estimated total error variance of 0.0038, which is a sum of the measurement and spatial representation error variances.

3 Results

3.1 Overview

Through the assimilation of ice-core data with a prior ensemble that is constant in time, we produce a spatially complete, mean-annual Greenland temperature and precipitation reanalysis at 50-year resolution (Figs. 4 and 5). Here we focus on results for the late glacial and the HTM.

Figure 4Spatial pattern of the reanalysis mean for temperature (a, c) and precipitation (b, d). Panels (a) and (b) are averaged over 1000 years around the peak warmth in the Holocene (5.5–4.5 ka), while panels (c) and (d) are averaged over 5000 years in the late glacial (20–15 ka). Anomalies and fractions are with respect to the mean of 1850–2000 CE. Points show ice-core locations used for each reanalysis with closed circles indicating δ18O records and open circles indicating accumulation records. Grey stars show the locations of the EGRIP ice-core site, Summit, and South Dome, which are referenced in Figs. 5 and 10.

Figure 5Time series of the prior (a, c, e) and reanalysis (b, d, f) ensemble mean and 5th to 95th percentile shading for temperature (red) and precipitation (blue) at three locations. Anomalies and fractions are with respect to the mean of 1850–2000 CE. Panels (a) and (b) show these time series for the location closest to the EGRIP ice-core site, panels (c) and (d) show them for the location closest to Summit, and panels (e) and (f) show for the location closest to South Dome. These locations are ordered from northernmost (top) to southernmost (bottom) and are shown on a map in Fig. 4.


In our reanalysis, late-glacial (20–15 ka) mean-temperature anomalies range from about −20C in northern Greenland to less than −10C in southern Greenland (Fig. 4c). At the GRIP and GISP2 ice-core sites, the reanalysis has a −14C anomaly with a standard deviation of 2 C. This is in excellent agreement with the mean-temperature anomaly of −14C for the same period at the GISP2 site, which was derived from δ18O calibrated with borehole thermometry (Cuffey et al.1995; Cuffey and Clow1997). Also in agreement with borehole thermometry, this period is warmer than the Last Glacial Maximum (Dahl-Jensen et al.1993; Cuffey and Clow1997). Average late-glacial precipitation in the reanalysis ranges from a third to half of modern with the highest values on the coasts around southern Greenland (Fig. 4d).

Our reanalysis shows warmest temperatures occurred across Greenland between 7 and 3 ka, reaching a temperature maximum around 5 ka (Fig. 5). Although this timing tends to be later than many estimates of the HTM, it lies within the ranges reported in the literature; for example, a summary of proxy records from around Greenland shows peak warmth usually occurring around 9–5 ka, but also as early as 10.8 ka and as late as 3 ka (Kaufman et al.2004). Borehole thermometry shows that temperatures peaked around 6–7.7 ka at Summit and 4.5 ka at the Dye3 ice-core site (Cuffey and Clow1997; Dahl-Jensen et al.1998; Kaufman et al.2004). Temperature estimates from δ15N of N2 show an earlier HTM peak at Summit around 8 ka (Kobashi et al.2017). In northwest Greenland, δ18O measurements from lake sediments show the HTM starting before 7.7 ka and ending around 6 ka (Lasher et al.2017), while chironomid assemblages show peak warmth around 10–8 ka (McFarlin et al.2018).

Mean-annual HTM temperature anomalies in our reanalysis range from nearly +2C in northern Greenland to about +1C in southern Greenland (Fig. 4a). Similar to the late-glacial temperature anomalies, the pattern of the HTM is dominated by a north–south trend that has the greater temperature changes to the north, especially in northwest Greenland. While this spatial pattern agrees well with previous studies which have noted especially warm Holocene temperature anomalies in northwest Greenland (Lasher et al.2017; Lecavalier et al.2017; McFarlin et al.2018), many estimates of HTM anomalies around Greenland are higher than our reanalysis indicates. Our low temperature estimates (compared to previous work) may be in part due to our reconstructing the annual mean rather than the summer mean; the greatest temperature anomalies in the HTM are thought to have occurred in the summer months, and many proxies are more sensitive to summer than annual temperature (Kaufman et al.2004). Importantly, in our reanalysis, the higher HTM temperatures do not translate into a marked increase in precipitation as would be expected from a thermodynamically scaled relationship between temperature and precipitation. Instead, we find fractional precipitation within 2 % of modern values (1.0±0.02) during the HTM (Fig. 4b), with slightly higher-than-modern precipitation in central Greenland and slightly lower-than-modern precipitation in northwestern Greenland.

3.2 Independent proxy evaluation

Here we evaluate our results against proxy records that are excluded from 10 of the iterations that make up the temperature and precipitation reanalyses. For this evaluation, we use the raw results (without a mean-bias correction). We find, however, that the mean biases are small relative to climate changes over the last 20 000 years; there is little difference between our bias-corrected and uncorrected results, and it is unlikely that the mean bias has a large effect on our evaluation.

Evaluation against independent proxy data shows that our reanalysis captures the timing and magnitude of low-frequency climate changes (Figs. 6 and 7) and is an improvement over both the prior ensemble and TraCE-21ka (Figs. S6–S7 and S8–S9). Evaluation over the full 20 000 years of the temperature and precipitation results shows high, positive correlation coefficients (ranging from a minimum of 0.97 to maximum of 0.99), which indicate that the reanalysis captures both the timing and sign of climate events, while high CE (0.87–0.98) and low RMSE values (0.62 ‰–1.2 ‰ for δ18O and 0.04–0.08 for accumulation) indicate that the reanalysis captures the magnitude of these events. The skill during this longest evaluation period is primarily due to the presence of low-frequency climate changes, which tend to be coherent across Greenland, such that evaluation over this full 20 000-year period shows more skill than evaluation over the full Holocene, which shows more skill than evaluation over just 5000 years in the Holocene (or 5000 years in the glacial) (Figs. 6 and 7).

Figure 6Skill metrics averaged over iterations and time for the temperature reanalysis. The first column (a, d, g, j) shows the skill metrics for the full overlap (Full) between the proxy record and reanalysis. A white dot indicates evaluation against proxy records that overlap only the Holocene (11.7–0 ka). The middle column (b, e, h, k) shows the skill metrics for a period in the glacial (Gl.) (20–15 ka), while the right column (c, f, i, l) is for a period in the Holocene (Hol.) (8–3 ka). The first row (a–c) reports the correlation coefficient, the second row (d–f) the coefficient of efficiency (CE), the third (g–i) the root mean square error (RMSE), and the fourth row (j–l) the ensemble calibration ratio (ECR). Triangle symbols pointing up indicate that the posterior ensemble evaluates better than the prior ensemble for that location and statistic. Triangle symbols pointing down indicate the opposite. We define better evaluation as correlation coefficient closer to 1, CE closer to 1, RMSE closer to 0, and ECR closer to 1.

Figure 7Skill metrics averaged over iterations and time for the precipitation reanalysis. The first column (a, d, g, j) shows the skill metrics for the full overlap (Full) between the proxy record and reanalysis. A white dot indicates evaluation against proxy records that overlap only the Holocene (11.7–0 ka). The middle column (b, e, h, k) shows the skill metrics for a period in the glacial (Gl.) (20–15 ka), while the right column (c, f, i, l) is for a period in the Holocene (Hol.) (8–3 ka). The first row (a–c) reports the correlation coefficient, the second row (d–f) the coefficient of efficiency (CE), the third (g–i) the root mean square error (RMSE), and the fourth row (j–l) the ensemble calibration ratio (ECR). Triangle symbols pointing up indicate that the posterior ensemble evaluates better than the prior ensemble for that location and statistic. Triangle symbols pointing down indicate the opposite. We define better evaluation as correlation coefficient closer to 1, CE closer to 1, RMSE closer to 0, and ECR closer to 1.

The posterior ensemble consistently shows improvement over the uninformed, constant prior ensemble during the 20 000-year evaluation period (Figs. S6 and S7). The TraCE-21ka simulation is also uninformed by the ice-core data, but it is transient and generally captures glacial to Holocene changes. Over the 20 000-year evaluation period, relative to the reconstructions, we find that TraCE-21ka has consistently lower correlation coefficients (0.86–0.96), lower CE values (0.50–0.86), and higher RMSE values (1.9 ‰–2.8 ‰ for δ18O and 0.11–0.15 for accumulation) (Figs. S8 and S9). This comparison suggests that our reconstruction captures the timing and magnitude of the glacial to Holocene transition better than TraCE-21ka.

Even for the shorter evaluation periods, which are dominated by high-frequency, spatially incoherent noise, the reanalysis shows improvement over both the prior ensemble and TraCE-21ka (Figs. S6–S7 and S8–S9); however, the improvement is not as consistent as for the 20 000-year evaluation period. For our reconstruction, correlation is positive except for three locations in the temperature reconstruction, with Holocene precipitation showing the largest correlation values (up to 0.60) and the total range being −0.30 to 0.60. The prior correlation is zero for these shorter evaluation periods and locations; however, TraCE-21ka shows correlation values ranging from −0.29 to 0.69, with eight negative correlations (more than we find for our reconstruction) but generally higher correlations in the Holocene than our reconstruction.

For the shorter evaluation periods, the reconstruction CE is generally negative (ranging from −82 to 0.17) with a few exceptions; however, the reconstruction may still be skillful (e.g., Cook et al.1999). The skill of the reconstruction is better measured by the difference between prior and posterior CE due to the strong influence that bias can have on CE (Hakim et al.2016). There is consistent improvement of the posterior CE over that of the prior (ranging from an increase of 4.7 to 3200) and over that of TraCE-21ka (ranging from a decrease of −6.9 to an increase of 230). RMSE is the most consistent of the skill metrics for these shorter evaluation periods, with our reconstruction showing improvement over both the prior and TraCE-21ka, the one exception being that TraCE-21ka has greater skill at the Agassiz ice-core site in the Holocene. For the reconstruction, RMSE values range from 0.24 ‰2 to 1.8 ‰2 for temperature and 0.025 to 0.10 for precipitation.

For all evaluation periods and both variables, the ECR for the prior is skewed towards values greater than unity (0.66–8.7), which suggests that the prior ensemble tends to have too little spread. Conversely, for the posterior, the ECR is generally less than unity (0.10–1.7) (Figs. 6 and 7), suggesting that the posterior ensemble has more spread than the error in the ensemble mean (as compared to the proxy records). This result indicates that the reanalysis ensemble encompasses the climate as recorded by the proxy records for most times and locations over the last 20 000 years.

3.3 Sensitivity evaluation

Proxy networks that change in time, such as ours, can introduce artificial discontinuities into the reanalysis, especially if the number of proxies is low or the proxy uncertainty values are inappropriate. We produce another reconstruction with a fixed proxy network, in which all assimilated proxy records participate in every time step in the reconstruction (see Sect. 2.3). A comparison of these results with our main reanalysis shows no apparent discontinuities for the ensemble mean and 5th to 95th percentiles at Summit (Fig. 8) or other locations around Greenland.

Figure 8Changing (red and blue) vs. fixed (black) proxy network for the (a) temperature (T) and (b) precipitation (P) reanalysis mean and 5th to 95th percentile shading. Anomalies and fractions are with respect to the mean of 1850–2000 CE. These time series are for the location closest to Summit, which is representative of the results around Greenland.


Our results are sensitive to the δ18OTsite relationship. To test this, we compare the main reanalysis to the four scenarios (S1–S4, as described in Sect. 2.3.1) that use different δ18O PSMs, each of which assumes a different slope and spatial pattern of the δ18OTsite relationship. We show this comparison for Summit as an example (Fig. 9), but the findings are applicable for all locations. As discussed previously, scenarios S1–S3 all assume that the δ18OTsite relationship has a uniform spatial pattern, but they each assume a different influence of precipitation seasonality. From these scenarios, we find that the temperature reconstruction is sensitive to the assumed precipitation seasonality, especially in the glacial where a stronger seasonality results in a greater glacial temperature anomaly. At Summit, this difference is nearly 10 C between S1, which assumes no influence, and S3, which assumes the most influence of precipitation seasonality (Fig. 9). Similarly, the main reanalysis and S4 scenario both assume that the δ18OTsite relationship has a spatially variable pattern, but S4 assumes a greater influence of precipitation seasonality. Again we find that the results are sensitive to assumed seasonality, with the greatest impact on the glacial-to-interglacial change.

Figure 9The main temperature (T) reanalysis (ensemble mean and 5th to 95th percentile shading) and ensemble mean for four sensitivity scenarios: S1–S4. Panel (a) shows the full 20 000-year reconstruction, while panel (b) shows a the Younger Dryas to early Holocene period (13 to 7 ka). Each sensitivity scenario reflects a different assumption about precipitation seasonality, with S1–S3 assuming a spatially uniform seasonality and S3–S4 assuming stronger seasonality than the main reanalysis. Anomalies are with respect to the mean of 1850–2000 CE. These time series are for the location closest to Summit, which is representative of the results around Greenland.


The temperature results are also sensitive to the spatial pattern of the δ18OTsite relationship. We find this by comparing the results from the S1–S3 scenarios that assume a spatially uniform relationship to results from the main reanalysis and S4 scenario that assume a spatially variable relationship. The S1–S3 scenarios have a characteristic shape to their time series (Fig. 9), and, although the main reanalysis and S4 scenario generally fit this characteristic shape in the glacial and middle–late Holocene, in the early Holocene the main reanalysis and S4 diverge and show slower warming trends than the S1–S3 scenarios (Fig. 9b). In addition, the reconstructions with spatially varying δ18OTsite relationships show stronger north–south gradients during times of abrupt temperature change than those with spatially constant relationships (Table S2). These findings indicate that there is new information added by using a PSM that includes spatial variability in precipitation seasonality.

For precipitation, we find that the results are sensitive to which accumulation record is assimilated at each ice-core site. As explained in Sects. 2.1 and S1, we use a low-, moderate-, and high-accumulation record for most of the ice-core sites to produce the low-, main-, and high-precipitation scenarios, respectively (Fig. 10). The largest spread among the scenarios is in the earliest part of the reconstruction, i.e., the last glacial through the early Holocene, since uncertainties in the ice-thinning history have the greatest impact at greater depths (and hence, greater ages). There is also a larger spread among the scenarios at more southern locations because the accumulation record at Dye3 has both the most influence on southern Greenland (Sect. S4) and the largest uncertainty in the ice-thinning history (Fig. 3, Sect. S2).

Figure 10Ensemble mean and 5th to 95th percentile shading for the main precipitation reanalysis (black), high-sensitivity scenario (blue), and low-sensitivity scenario (red). Fractions are with respect to the mean of 1850–2000 CE. Panel (a) is the time series for the location closest to the EGRIP ice-core site, panel (b) is closest to Summit, and panel (c) is closest to South Dome, which are representative of northern, central, and southern Greenland and are shown on a map in Fig. 4.


4 Discussion

4.1 Precipitation–temperature relationship

Our results allow us to investigate the relationship between temperature and precipitation. To facilitate comparison with thermodynamic scalings widely used by ice-sheet models, we define the relationship as exponential and find the best fit for our reanalysis:

(8) P fraction = P past P modern = e β ( T past - T modern ) ,

where P is the precipitation rate, T is the temperature, and β is a scaling factor (Greve et al.2011). For a given temperature change, a higher value of β results in a larger change in precipitation (orange in Fig. 11a). In ice-sheet models that use this scaling, it is commonly applied with a uniform β value for all locations (e.g., Huybrechts et al.1991; Huybrechts2002; Greve1997; Greve et al.2011; Pollard and DeConto2012; Cuzzone et al.2019). The best-fit scaling factors (β) center on the Greve et al. (2011) value of 0.07 for locations around Greenland, but the scaling factors tend to be larger (β>0.10) where late-glacial precipitation is lowest and smaller (β<0.05) where late-glacial precipitation is highest (Figs. 4d and 11a).

Figure 11The precipitation–temperature relationship in our main reanalysis. Panel (a) shows the spatial pattern of the scaling factor (β) for the best-fit thermodynamic scaling. The color bar is centered on 0.07, the value used by Greve et al. (2011). Points indicate ice-core locations used for each reanalysis with closed circles indicating δ18O records and open circles indicating accumulation records. The star is at the center of the area used in panels (b)(e) (65 to 68.7° N and 48.5 to 52.5° W). Panel (b) is a scatterplot of temperature anomaly vs. precipitation fraction from the reanalysis (blue points). The black line shows the best-fit exponential scaling. Panel (c) shows the time series of the precipitation reanalysis (blue line) and precipitation scaled from temperature using the best-fit scaling (black line). Panels (d) and (e) are the same as panel (c) except low-pass- and high-pass-filtered, respectively, with a cutoff frequency of 5000 per year.

Spatial variability in the temperature–precipitation relationship was also found by a previous study, which looked at decadal averages over a recent 110-year period (Buchardt et al.2012). Also using an exponential to represent the relationship, Buchardt et al. (2012) similarly found weaker relationships in northern than in central and southern Greenland (compare their Fig. 3 to our Fig. 11). In southern Greenland, they found a negative β value in southwest and a positive value in southeast Greenland, which they attribute to the Foehn effect. We find the opposite pattern, with larger, positive β values in southwest Greenland as compared to southeast Greenland. These differences may result from the time periods, timescales, and spatial resolutions considered in each investigation and from methodology. The spatial pattern found by Buchardt et al. (2012) is limited to the ice-core locations they considered such that they have little information about the relationship near the ice-sheet edges. The spatial pattern in our reanalysis is influenced by both the ice-core records and the spatial covariance structures of the prior ensembles. With only one long ice core in southern Greenland, the east–west gradient in β values is highly influenced by the spatial patten inherited from TraCE-21ka.

The precipitation–temperature relationship in our reanalysis is driven by the assimilated ice-core records, though, as mentioned above, the spatial pattern of this relationship is also influenced by the spatial covariance structures of the prior ensembles. Previous work with ice-core records has found that the relationship between temperature and precipitation is frequency dependent, with a stronger relationship at lower frequencies (Cuffey and Clow1997); as expected, there is a similar frequency dependence in our reanalysis. We find that an exponential scaling captures the low-frequency glacial to Holocene precipitation change; however, the same scaling fails at higher frequencies (Fig. 11b–e). To evaluate this frequency dependence, we filtered our results using sixth-order, low-pass and high-pass Butterworth filters with 5000 per year cutoff frequencies. The low-pass-filtered dataset shows the same precipitation–temperature relationship as the unfiltered dataset (Fig. 11c–d), while the high-pass-filtered dataset shows that precipitation is less sensitive to changes in temperature (i.e., the value of β is lower) at these higher frequencies (Fig. 11e). This decoupling of temperature and precipitation is apparent in the amplitude difference of high-frequency signals in the glacial and the Holocene. In our temperature reanalysis, we find that high frequencies in the glacial have a greater amplitude than those in the Holocene, while in our precipitation reanalysis, we find the opposite. A single scaling, as is typically used in ice-sheet modeling, cannot capture this difference.

We examine how the sensitivity experiments (Figs. 9 and 10) affect the scaling factor (β) in the precipitation–temperature relationship. We pair the five temperature reconstructions (main, S1–S4) and three precipitation reconstructions (low, moderate, and high) into 15 possible combinations and conduct the same analysis as described above. Across these 15 combinations, we find that the spatial pattern of β is robust (Fig. 11a). The exact magnitude depends primarily on the temperature reconstruction and how cold it is in the glacial, with colder temperatures giving lower β values. To a lesser degree, the magnitude also depends on the precipitation reconstruction, with wetter scenarios giving lower β values. As previously, we find that the low-pass-filtered datasets have the same or nearly the same β value as the unfiltered dataset, while the high-pass-filtered datasets have consistently lower β values. As an example of this, Table S3 shows the β value found for the Kangerlussuaq region for all 15 combinations and three filtering options.

4.2 Spatial patterns during abrupt climate change events

Our paleoclimate data assimilation framework is not limited to the assimilation of δ18O and accumulation rate records. In this section we examine how our reanalysis compares to reconstructions driven by another type of proxy: δ15N of N2. In particular, we focus on abrupt temperature events, for which there is previous work using δ15N (e.g., Severinghaus et al.1998; Guillevic et al.2013; Buizert et al.2014). Abrupt temperature events increase the thermal gradient in the firn – the upper, porous portion of the ice column – which leads to fractionation of the stable isotopes of N2 (Severinghaus et al.1998). Using a firn compaction model, temperature can be derived from δ15N measurements with inverse methods (e.g., Severinghaus et al.1998; Severinghaus and Brook1999; Guillevic et al.2013; Kindler et al.2014; Buizert et al.2014; Kobashi et al.2017). We assimilate temperatures derived from δ15N data collected from the GISP2, NGRIP, and NEEM ice cores (Buizert et al.2014).

Our reanalysis and the Buizert et al. (2014) records cover three abrupt temperature events – the Bølling–Allerød warming (14.6 ka), cooling into the Younger Dryas (12.9 ka), and warming at the end of the Younger Dryas (11.5 ka). For these three abrupt temperature events, the δ15N-derived temperature records show larger temperature changes at GISP2 in central Greenland (Buizert et al.2014), while the δ18O changes are larger at NEEM in northwest Greenland (Fig. 2). We perform three sets of experiments to investigate how these two sets of proxy records affect the mean spatial pattern indicated by our reanalysis during abrupt temperature events. The first experiment, “O8”, assimilates all eight δ18O records with one 100-member prior ensemble from TraCE-21ka. This results in a 100-member reanalysis ensemble. The second experiment, “N3O5”, assimilates all three of the δ15N-derived temperature records and the five remaining δ18O records (those that do not overlap with the δ15N sites), using the same 100-member prior ensemble as used in the O8 experiment. Finally, we perform a modified experiment, “N3O5_BA”, with both δ18O and δ15N records, but using a prior ensemble selected from the 1000 years surrounding the Bølling–Allerød warming. To maintain a 100-member prior ensemble, we use decadal rather than 50-year averages of TraCE-21ka for this experiment. This adjustment does not affect the comparison (results not shown). Detailed methods for these experiments are given in Sect. S6.

For both the “O8” and “N3O5” experiments, the spatial pattern for the abrupt climate change events are similar to our main reanalysis, with the largest magnitude temperature changes in northern and northwestern Greenland, decreasing magnitude with decreasing latitude, and slightly larger change in the central east and southeast than corresponding western regions. For example, the spatial pattern of the Younger Dryas cooling is nearly the same regardless of which grouping of records is assimilated (Fig. 12a and b). This overall finding is robust to different combinations of these proxy records; for example, if we assimilate just the three δ15N-derived temperature records and no δ18O records (results not shown) the pattern is not substantially changed. This pattern of temperature change differs from spatial patterns inferred previously from δ15N for various abrupt climate events (Guillevic et al.2013; Buizert et al.2014); however, the O8 and N3O5 experiments show that these differences are not due to the assimilation of δ18O rather than δ15N-derived temperatures. The replacement of δ18O records with δ15N-derived temperature records does not change the overall spatial pattern of these abrupt events, though it does result in a reconstruction with slightly larger temperature changes in the south and slightly smaller changes in the north. This effect, however, is less important than the choice of prior ensemble, as is shown by the third experiment, N3O5_BA.

Figure 12Spatial pattern of the abrupt cooling event into the Younger Dryas. Panel (a) shows results from experiment O8, assimilating all eight δ18O records, panel (b) shows results from experiment N3O5, assimilating all three δ15N-derived temperature records and the remaining five δ18O records (those that do not overlap with the δ15N sites), and panel (c) shows results from experiment N3O5_BA, which is similar to the N3O5 experiment except that the prior ensemble is selected from the 1000 years surrounding the Bølling–Allerød warming. Unfilled black circles show locations of assimilated δ18O records, while filled circles with white outlines show locations of assimilated δ15N-derived temperature records. Filled circles in panels (b) and (c) show the δ15N-derived temperature values as reported by Buizert et al. (2014) on the same color scale as the rest of the panel. The temporal definition of this event is the same as defined in Buizert et al. (2014).

For N3O5_BA, we restrict the prior ensemble to the Bølling–Allerød warming, which produces a reconstructed spatial pattern similar to those reported by Buizert et al. (2014) for each of the three abrupt climate events. In TraCE-21ka, the Bølling–Allerød warming is forced by a sudden termination of freshwater forcing in the North Atlantic (Liu et al.2009). This forcing leads to large temperature fluctuations in southern Greenland that decrease with increasing latitude. With this covariance pattern dominating the prior ensemble, the N3O5_BA reconstruction does indeed show the largest temperature changes in southern Greenland, followed by central and then northern Greenland (Figs. 12c, S12c, and S13c).

These experiments suggest that current ice-core records are insufficient to place a strong constraint on the spatial pattern of abrupt climate events. Reconstructions of these events may be improved with additional data, especially from southern Greenland (Sect. S4), and with prior ensembles that are designed to sample over uncertainty in the forcing and boundary conditions unique to these events.

4.3 Climate and the ice sheet: a case study of southwest Greenland

Our main reanalysis is one of very few spatially complete time series of Greenland climate over the last 20 000 years. Here, we compare our reanalysis with other Greenland climate histories, and suggest that, together, they should be treated as an ensemble of climate boundary conditions that can be used to produce ensembles of ice-sheet model simulations. These climate histories can also be further evaluated using a combination ice-sheet models and independent constraints from the glacial–geologic record of past ice-sheet configurations.

We compare our results with the recent reconstruction of Buizert et al. (2018), hereafter referred to as the B18 reconstruction. The B18 temperature reconstruction was produced by adjusting a part of the TraCE-21ka temperature field that is affected by changes in the Atlantic Meridional Overturning Circulation such that the full temperature field provides a good match to an average of three δ15N-derived temperature records recovered from ice cores. The B18 snow-accumulation reconstruction is simply a reference climatology scaled to accumulation rates from the GISP2 ice core. We treat this as a precipitation reconstruction but note that accumulation may be less than precipitation at some locations around Greenland, especially near the coast. It is also informative to compare these results with our S4 temperature and high-precipitation reconstructions (hereafter referred to as S4 and high P or simply “sensitivity” as in Fig. 13), as well as the TraCE-21ka simulation itself (i.e., the climate-model output, unconstrained by data). For brevity, we focus on the area around Kangerlussuaq in southwest Greenland, but the comparisons are generally applicable to any region of Greenland. Southwest Greenland is of interest because the ice-sheet behavior here is primarily a response to changes in surface forcing (i.e., temperature and precipitation) because there are few tidewater glaciers (Cuzzone et al.2019). Furthermore, the Kangerlussuaq region has a particularly well-documented ice-sheet retreat history through the Holocene (Young and Briner2015; Lesnek and Briner2018).

Figure 13Temperature (a) and precipitation (b) reconstructions from our main reanalysis (black), our sensitivity scenarios S4 and high P (green), TraCE-21ka (red), and B18 (blue) (Buizert et al.2018). Each reconstruction is averaged to a 50-year time resolution and averaged over a spatial domain in the Kangerlussuaq region, defined by the latitude-longitude box 65 to 68.7° N and 48.5 to 52.5° W, the center of which is located at the star in Fig. 11a. Temperature anomalies and precipitation fraction are defined with reference to the mean of 1850–2000 CE.


In the Kangerlussuaq region, the B18 reconstruction shows more extreme temperature changes than our reconstructions, with late-glacial (20–15 ka) anomalies of about −20C and peak HTM temperature anomalies of about +2C at 9 ka (Fig. 13). B18 also shows a faster rate of transition between the glacial and Holocene, reaching temperatures close to modern values by 10 ka. In contrast, TraCE-21ka shows more moderate temperature anomalies and a slower transition, with late-glacial anomalies of about −8.6C and near-modern temperatures that first appear around 7 ka. TraCE-21ka has no obvious HTM in this location or any location in Greenland. Our main reanalysis and the S4 version of our temperature reanalysis both lie between B18 and TraCE-21ka, with late-glacial anomalies of about −12 and −14C, respectively, Holocene peak temperature anomalies of +1C around 5 ka, and temperatures close to modern values first appearing around 8 ka.

For precipitation, the B18 reconstruction again tends to show the largest fluctuations and fastest transition, with a late-glacial precipitation fraction of about 0.26 and precipitation rates close to modern values first appearing just after 10 ka. TraCE-21ka again shows the most moderate fluctuations and a slower transition, with a late-glacial fraction of about 0.38 and rates close to modern values not appearing until around 5 ka. Our main reanalysis and high P lie in the middle during the late glacial, with fractions of about 0.32 and 0.36, respectively; however, our main reanalysis has a slow transition into the Holocene, similar to TraCE-21ka, while high P has a fast transition similar to B18. In the Holocene, high P shows the most elevated precipitation out of all the reconstructions, with 10 %–15 % more precipitation than modern occurring around 7–3 ka. B18 shows precipitation values similar to modern values for the last 10 000 years of the Holocene, while TraCE-21ka and our main reanalysis show lower-than-modern precipitation throughout most of the Holocene.

All of these paleoclimate reconstructions – our main reanalysis, the sensitivity scenarios, and B18 – are plausible histories of temperature and precipitation over Greenland. Given any past change in the ice sheet, each of these histories has a different implication for ice-sheet sensitivity to climate, the veracity of which could be tested by using them to force an ice-sheet model and comparing this ensemble of results to the geologic record.

Our results have potentially important implications for the response of the Greenland Ice Sheet to climate change. In particular, we find maximum Holocene temperatures were reached around 5 ka, which is between 500 years earlier and 4000 years later than most previous estimates. Moreover, there is little corresponding change in precipitation in our main reanalysis and the low-sensitivity scenario. If these findings are correct, they imply a relatively rapid response to temperature forcing for sections of the ice sheet margin that retreated less than a century later (Young and Briner2015). A caveat is that proxy data remain very sparse, particularly in southern Greenland, where the poorly resolved Dye3 core is the only long record. Future work to obtain improved measurements on the Dye3 core or to gather new data from southern Greenland would help to alleviate this limitation, as would the incorporation of data from off the ice sheet, such as from lake and ocean sediment cores.

An important distinction among various different paleoclimate reconstructions for Greenland is in the treatment of elevation changes. Any paleoclimate reconstruction from ice-core records is complicated by ice-sheet elevation changes. In Vinther et al. (2009), it is assumed that the climate history is the same at all locations around Greenland and that any differences among the ice-core paleotemperature records is a result of that elevation change. In B18, past elevation changes are assumed to be negligible. In our reconstruction, the impact of elevation change on the spatial covariances of temperature and precipitation is implicitly accounted for as part of the data assimilation methodology. Formally, our reconstruction is of surface climate, not climate at a fixed elevation. Consequently, our reanalysis may not be directly comparable to other paleoclimate reconstructions. For example, the HTM is commonly reconstructed as an early Holocene event in records that are at a fixed or nearly fixed elevation. In our reanalysis, the early Holocene is cooler than the mid Holocene. Changes in the ice-surface elevation could account for this apparent discrepancy. Thinning in the early Holocene (Vinther et al.2009) would result in a lowering of the ice surface and an apparent warming at the ice surface due to lapse rate effects. This warming signal would be captured in ice-core records. If the warming trend due to surface lowering occurs at the same time as an overall climate cooling, then the climate signal would be dampened or possibly reversed.

Our method depends on the accuracy of the climate–elevation relationships in our prior – i.e., in the TraCE-21ka climate-model simulation, which probably does not capture such relationships with particularly high fidelity since the model resolution is low and the climate and ice-sheet models are not coupled. Future work could take advantage of the probabilistic relationships among accumulation, temperature, and surface elevation as simulated in fine-scale regional climate models (Edwards et al.2014).

5 Conclusions

Paleoclimate data assimilation is a novel method for reconstructing climate fields over the Greenland Ice Sheet. Our approach, combining ice-core records with a climate-model simulation, provides complete spatial reanalyses of mean-annual temperature and precipitation covering the last 20 000 years. Evaluation against independent proxy records shows that this methodology leads to significant and meaningful improvement over the prior ensemble (drawn from a climate simulation) and TraCE-21ka. Between the posterior ensemble and sensitivity experiments, our results provide a range of climate scenarios for ice-sheet modeling. Moreover, independently reconstructing both precipitation and temperature allows the assumption of purely thermodynamic control on precipitation to be relaxed and an examination of the relationship between these quantities over a range of timescales. Specifically, we find that the Clausius–Clapeyron scaling is a good approximation over glacial–interglacial cycles but not for shorter timescales, where precipitation variability partially decouples from temperature.

Paleoclimate reconstructions would benefit from a larger selection of long climate-model simulations at higher resolution. Particularly valuable would be transient simulations that include water isotopes as prognostic variables, which allow for direct assimilation of water isotope ratios (Steiger et al.2017; Okazaki and Yoshimura2017), rather than the use of an explicit proxy system model between temperature and δ18O. Recent work shows significant improvements to the realism of water-isotope enabled models in the polar regions (Nusbaumer et al.2017; Dütsch et al.2019; Cauquoin et al.2019; Okazaki and Yoshimura2019), and longer simulations, once available, should allow us to further improve upon the results we have presented here. In principle, our method could also be applied to climate-model simulations that include a fully coupled Greenland Ice Sheet. At present, fully coupled simulations of Greenland over thousands of years are prohibitively expensive except at low resolution, and the limited work that has been done shows significant biases (Vizcaino et al.2015). Nevertheless, incorporating data assimilation into such models would provide the groundwork for more-complete data-constrained simulations as computing power becomes less of a limiting factor in the future.

Code and data availability

The paleoclimate reconstructions in this paper made use of code from the Last Millennium Reanalysis project, which is publicly available at (Hakim2019). The reconstructions and ice-core records, including the new accumulation histories for Dye3, GRIP, and NGRIP, are publicly available at the Arctic Data Center (, Badgeley et al.2020).


The supplement related to this article is available online at:

Author contributions

JAB and EJS conceived the idea for the study. JAB wrote code improvements necessary for this work, completed the paleoclimate reconstructions with guidance from EJS and GJH, conducted all analyses of the results, and wrote the first draft of the paper. GJH provided expert advice on data assimilation methodology and code development. TJF made the calculations of ice flow used to model accumulation for the Dye3, GRIP, and NGRIP cores. All authors contributed to the final version of the paper.

Competing interests

The authors declare that they have no conflict of interest.


Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation.


We thank Robert Tardif for help with code development, and we thank Joshua Anderson, Bo Vinther, and Christo Buizert for their help compiling the ice-core data. We also thank the Snow on Ice project members for their discussions and support, especially Joshua Cuzzone, Jason Briner, and Elizabeth Thomas.

Financial support

This research has been supported by the National Science Foundation, Office of Polar Programs (grant no. 1503281), the National Science Foundation, Division of Graduate Education (grant no. 1256082), the Heising Simons Foundation (grant no. 2016-14), and the National Science Foundation, Division of Atmospheric and Geospace Sciences (grant no. 1602223).

Review statement

This paper was edited by Amaelle Landais and reviewed by three anonymous referees.


Alley, R. B., Shuman, C., Meese, D., Gow, A., Taylor, K., Cuffey, K., Fitzpatrick, J., Grootes, P., Zielinski, G., Ram, M., Spinelli, G., and Elder, B.: Visual-stratigraphic dating of the GISP2 ice core: Basis, reproducibility, and application, J. Geophys. Res.-Oceans, 102, 26367–26381, 1997. a

Alley, R. B., Andrews, J. T., Brigham-Grette, J., Clarke, G., Cuffey, K. M., Fitzpatrick, J., Funder, S., Marshall, S., Miller, G., Mitrovica, J., Muhs, D., Otto-Bliesner, B., Polyak, L., and White, J.:: History of the Greenland Ice Sheet: paleoclimatic insights, Quaternary Sci. Rev., 29, 1728–1756, 2010. a

Andersen, K. K., Azuma, N., Barnola, J.-M., Bigler, M., Biscaye, P., Caillon, N., Chappellaz, J., Clausen, H. B., Dahl-Jensen, D., Fischer, H., Flückiger, J., Fritzsche, D., Fujii, Y., Goto-Azuma, K., Grønvold, K., Gundestrup, N. S., Hansson, M., Huber, C., Hvidberg, C. S., Johnsen, S. J., Jonsell, U., Jouzel, J., Kipfstuhl, S., Landais, A., Leuenberger, M., Lorrain, R., Masson-Delmotte, V., Miller, H., Motoyama, H., Narita, H., Popp, T., Rasmussen, S. O., Raynaud, D., Röthlisberger, R., Ruth, U., Samyn, D., Schwander, J., Shoji, H., Siggard-Andersen, M. L., Steffensen, J. P., Stocker, T., Sveinbjörnsdottir, A. E., Svensson, A., Takata, M., Tison, J. L., Thorsteinsson, T., Watanabe, O., Wilhelms, F., and White, J.: High-resolution record of Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, 2004. a

Andersen, K. K., Svensson, A., Johnsen, S. J., Rasmussen, S. O., Bigler, M., Röthlisberger, R., Ruth, U., Siggaard-Andersen, M.-L., Steffensen, J. P., Dahl-Jensen, D., Vinther, B. M., and Clausen, H. B.: The Greenland ice core chronology 2005, 15–42 ka. Part 1: constructing the time scale, Quaternary Sci. Rev., 25, 3246–3257, 2006. a, b

Armengaud, A., Koster, R. D., Jouzel, J., and Ciais, P.: Deuterium excess in Greenland snow: Analysis with simple and complex models, J. Geophys. Res.-Atmos., 103, 8947–8953, 1998. a

Badgeley, J. A., Steig, E. J., Hakim, G. J., and Fudge, T. J.: Reconstructions of mean-annual Greenland temperature and precipitation for the past 20 000 years and the ice-core records used to create the reconstructions, Arctic Data Center,, 2020. a

Bindschadler, R. A., Nowicki, S., Abe-Ouchi, A., Aschwanden, A., Choi, H., Fastook, J., Granzow, G., Greve, R., Gutowski, G., Herzfeld, U., Jackson, C., Johnson, J., Khroulev, C., Levermann, A., Lipscomb, W., Martin, M., Morlighem, M., Parizek, B., Pollard, D., Price, S., Ren, D., Saito, F., Sato, T., Seddik, H., Seroussi, H., Takahashi, K., Walker, R., and Wang, W.: Ice-sheet model sensitivities to environmental forcing and their use in projecting future sea level (the SeaRISE project), J. Glaciol., 59, 195–224, 2013. a

Bintanja, R., van de Wal, R. S., and Oerlemans, J.: Modelled atmospheric temperatures and global sea levels over the past million years, Nature, 437, 125–128, 2005. a

Buchardt, S. L., Clausen, H. B., Vinther, B. M., and Dahl-Jensen, D.: Investigating the past and recent δ18O-accumulation relationship seen in Greenland ice cores, Clim. Past, 8, 2053–2059,, 2012. a, b, c

Buizert, C., Gkinis, V., Severinghaus, J. P., He, F., Lecavalier, B. S., Kindler, P., Leuenberger, M., Carlson, A. E., Vinther, B., MassonDelmotte, V., White, J., Liu, Z., Otto-Bliesner, B., and Brook, E.: Greenland temperature response to climate forcing during the last deglaciation, Science, 345, 1177–1180, 2014. a, b, c, d, e, f, g, h, i, j

Buizert, C., Keisling, B., Box, J., He, F., Carlson, A., Sinclair, G., and DeConto, R.: Greenland-Wide Seasonal Temperatures During the Last Deglaciation, Geophys. Res. Lett., 45, 1905–1914, 2018. a, b, c, d

Cauquoin, A., Werner, M., and Lohmann, G.: Water isotopes – climate relationships for the mid-Holocene and preindustrial period simulated with an isotope-enabled version of MPI-ESM, Clim. Past, 15, 1913–1937,, 2019. a, b

Conway, H., Hall, B., Denton, G., Gades, A., and Waddington, E.: Past and future grounding-line retreat of the West Antarctic Ice Sheet, Science, 286, 280–283, 1999. a

Cook, E. R., Meko, D. M., Stahle, D. W., and Cleaveland, M. K.: Drought reconstructions for the continental United States, J. Climate, 12, 1145–1162, 1999. a

Cuffey, K. M. and Clow, G. D.: Temperature, accumulation, and ice sheet elevation in central Greenland through the last deglacial transition, J. Geophys. Res.-Oceans, 102, 26383–26396, 1997. a, b, c, d, e, f, g, h, i, j, k, l

Cuffey, K. M. and Steig, E. J.: Isotopic diffusion in polar firn: implications for interpretation of seasonal climate parameters in ice-core records, with emphasis on central Greenland, J. Glaciol., 44, 273–284, 1998. a

Cuffey, K. M., Clow, G. D., Alley, R. B., Stuiver, M., Waddington, E. D., and Saltus, R. W.: Large arctic temperature change at the Wisconsin-Holocene glacial transition, Science, 270, 455–458, 1995. a, b, c

Cuzzone, J. K., Schlegel, N.-J., Morlighem, M., Larour, E., Briner, J. P., Seroussi, H., and Caron, L.: The impact of model resolution on the simulated Holocene retreat of the southwestern Greenland ice sheet using the Ice Sheet System Model (ISSM), The Cryosphere, 13, 879–893,, 2019. a, b

Dahl-Jensen, D., Johnsen, S., Hammer, C., Clausen, H., and Jouzel, J.: Past accumulation rates derived from observed annual layers in the GRIP ice core from Summit, Central Greenland, in: Ice in the Climate System, edited by: Peltier, W. R. , NATO ASI Series (Series I: Global Environmental Change), vol. 12, Springer, Berlin, Heidelberg, Germany, 1993. a, b, c

Dahl-Jensen, D., Mosegaard, K., Gundestrup, N., Clow, G. D., Johnsen, S. J., Hansen, A. W., and Balling, N.: Past temperatures directly from the Greenland ice sheet, Science, 282, 268–271, 1998. a, b, c, d, e

Dahl-Jensen, D., Gundestrup, N., Gogineni, S. P., and Miller, H.: Basal melt at NorthGRIP modeled from borehole, ice-core and radio-echo sounder observations, Ann. Glaciol., 37, 207–212, 2003. a, b, c

Dahl-Jensen, D., Albert, M., Aldahan, A., Azuma, N., Balslev-Clausen, D., Baumgartner, M., Berggren, A.-M., Bigler, M., Binder, T., Blunier, T., et al.: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494, 2013. a, b, c

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468, 1964. a

Dansgaard, W. and Johnsen, S.: A flow model and a time scale for the ice core from Camp Century, Greenland, J. Glaciol., 8, 215–223, 1969. a

Dansgaard, W., Clausen, H., Gundestrup, N., Hammer, C., Johnsen, S., Kristinsdottir, P., and Reeh, N.: A new Greenland deep ice core, Science, 218, 1273–1277, 1982. a

Dütsch, M., Blossey, P. N., Steig, E. J., and Nusbaumer, J. M.: Non-equilibrium fractionation during ice cloud formation in iCAM5: evaluating the common parameterization of supersaturation as a linear function of temperature, J. Adv. Model. Earth Syst., 11, 3777–3793, 2019. a

Edwards, T. L., Fettweis, X., Gagliardini, O., Gillet-Chaulet, F., Goelzer, H., Gregory, J. M., Hoffman, M., Huybrechts, P., Payne, A. J., Perego, M., Price, S., Quiquet, A., and Ritz, C.: Probabilistic parameterisation of the surface mass balance–elevation feedback in regional climate model simulations of the Greenland ice sheet, The Cryosphere, 8, 181–194,, 2014. a

Fisher, D. A., Reeh, N., and Clausen, H.: Stratigraphic noise in time series derived from ice cores, Ann. Glaciol., 7, 76–83, 1985. a, b

Fudge, T., Markle, B. R., Cuffey, K. M., Buizert, C., Taylor, K. C., Steig, E. J., Waddington, E. D., Conway, H., and Koutnik, M.: Variable relationship between accumulation and temperature in West Antarctica for the past 31 000 years, Geophys. Res. Lett., 43, 3795–3803, 2016. a

Gierz, P., Werner, M., and Lohmann, G.: Simulating climate and stable water isotopes during the L ast I nterglacial using a coupled climate-isotope model, J. Adv. Model. Earth Syst., 9, 2027–2045, 2017. a

Gkinis, V., Simonsen, S. B., Buchardt, S. L., White, J., and Vinther, B. M.: Water isotope diffusion rates from the NorthGRIP ice core for the last 16 000 years – Glaciological and paleoclimatic implications, Earth Planet. Sc. Lett., 405, 132–141, 2014. a, b

Greve, R.: Application of a polythermal three-dimensional ice sheet model to the Greenland ice sheet: response to steady-state and transient climate scenarios, J. Climate, 10, 901–918, 1997. a, b

Greve, R., Saito, F., and Abe-Ouchi, A.: Initial results of the SeaRISE numerical experiments with the models SICOPOLIS and IcIES for the Greenland ice sheet, Ann. Glaciol., 52, 23–30, 2011. a, b, c, d, e, f

Grootes, P. and Stuiver, M.: Oxygen 18/16 variability in Greenland snow and ice with 10- 3-to 105-year time resolution, J. Geophys. Res.-Oceans, 102, 26455–26470, 1997. a

Guillevic, M., Bazin, L., Landais, A., Kindler, P., Orsi, A., Masson-Delmotte, V., Blunier, T., Buchardt, S. L., Capron, E., Leuenberger, M., Martinerie, P., Prié, F., and Vinther, B. M.: Spatial gradients of temperature, accumulation and δ18O-ice in Greenland over a series of Dansgaard–Oeschger events, Clim. Past, 9, 1029–1051,, 2013. a, b, c, d

Hakim, G. J.: Source code for the Last Millennium Reanalysis (LMR) project, production_develop, 7e6da55, available at:, last access: 25 April 2019. a

Hakim, G. J., Emile-Geay, J., Steig, E. J., Noone, D., Anderson, D. M., Tardif, R., Steiger, N., and Perkins, W. A.: The last millennium climate reanalysis project: Framework and first results, J. Geophys. Res.-Atmos., 121, 6745–6764, 2016. a, b, c, d

Harrison, S., Bartlein, P., Brewer, S., Prentice, I., Boyd, M., Hessler, I., Holmgren, K., Izumi, K., and Willis, K.: Climate model benchmarking with glacial and mid-Holocene climates, Clim. Dynam., 43, 671–688, 2014. a

Hawkins, E. and Sutton, R.: The potential to narrow uncertainty in regional climate predictions, B. Am. Meteorol. Soc., 90, 1095–1108, 2009. a

He, F., Shakun, J. D., Clark, P. U., Carlson, A. E., Liu, Z., Otto-Bliesner, B. L., and Kutzbach, J. E.: Northern Hemisphere forcing of Southern Hemisphere climate during the last deglaciation, Nature, 494, 81–85, 2013. a, b, c

Houtekamer, P. L., Mitchell, H. L., Pellerin, G., Buehner, M., Charron, M., Spacek, L., and Hansen, B.: Atmospheric data assimilation with an ensemble Kalman filter: Results with real observations, Mon. Weather Rev., 133, 604–620, 2005. a

Huybers, P. and Wunsch, C.: A depth-derived Pleistocene age model: Uncertainty estimates, sedimentation variability, and nonlinear climate change, Paleoceanography, 19, PA1028,, 2004. a

Huybrechts, P.: Sea-level changes at the LGM from ice-dynamic reconstructions of the Greenland and Antarctic ice sheets during the glacial cycles, Quaternary Sci. Rev., 21, 203–231, 2002. a, b, c

Huybrechts, P., Letreguilly, A., and Reeh, N.: The Greenland ice sheet and greenhouse warming, Palaeogeogr. Palaeocl., 89, 399–412, 1991. a, b

Johnsen, S., Dansgaard, W., and White, J.: The origin of Arctic precipitation under present and glacial conditions, Tellus B, 41, 452–468, 1989. a, b

Johnsen, S. J., Clausen, H. B., Dansgaard, W., Gundestrup, N. S., Hammer, C. U., Andersen, U., Andersen, K. K., Hvidberg, C. S., DahlJensen, D., Steffensen, J. P., Shoji, H., Sveinbjörnsdóttir, A. E., White, J., Jouzel, J., and Fisher, D.: The δ18O record along the Greenland Ice Core Project deep ice core and the problem of possible Eemian climatic instability, J. Geophys. Res.-Oceans, 102, 26397–26410, 1997. a

Jouzel, J., Alley, R. B., Cuffey, K., Dansgaard, W., Grootes, P., Hoffmann, G., Johnsen, S. J., Koster, R., Peel, D., Shuman, C., Stievenard, M., Stuiver, M., and White, J.: Validity of the temperature reconstruction from water isotopes in ice cores, J. Geophys. Res.-Oceans, 102, 26471–26487, 1997. a, b, c, d, e

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., and Joseph, D.: The NCEP/NCAR 40-year reanalysis project, B. Am. Meteorol. Soc., 77, 437–472, 1996. a

Kapsner, W., Alley, R. B., Shuman, C., Anandakrishnan, S., and Grootes, P.: Dominant influence of atmospheric circulation on snow accumulation in Greenland over the past 18 000 years, Nature, 373, 52–54, 1995. a

Kaufman, D., Ager, T., Anderson, N., Anderson, P., Andrews, J., Bartlein, P., Brubaker, L., Coats, L., Cwynar, L., Duvall, M., Dyke, A., Edwards, M., Eisner, W., Gajewski, K., Geirsdóttir, A., Hu, F., Jennings, A., Kaplan, M., Kerwin, M., Lozhkin, A., MacDonald, G., Miller, G., Mock, C., Oswald, W., Otto-Bliesner, B., Porinchu, D., Rühland, K., Smol, J., Steig, E., and Wolfe, B.: Holocene thermal maximum in the western Arctic (0–180 W), Quaternary Sci. Rev., 23, 529–560, 2004. a, b, c, d

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. a, b

Kobashi, T., Menviel, L., Jeltsch-Thömmes, A., Vinther, B. M., Box, J. E., Muscheler, R., Nakaegawa, T., Pfister, P. L., Döring, M., Leuenberger, M., Wanner, H., and Ohmura, A.: Volcanic influence on centennial to millennial Holocene Greenland temperature change, Sci. Rep., 7, 1441,, 2017. a, b, c

Krinner, G. and Werner, M.: Impact of precipitation seasonality changes on isotopic signals in polar ice cores: a multi-model analysis, Earth Planet. Sc. Lett., 216, 525–538, 2003. a

Langen, P., Mottram, R., Christensen, J., Boberg, F., Rodehacke, C., Stendel, M., Van As, D., Ahlstrøm, A., Mortensen, J., Rysgaard, S., Petersen, D., Svendsen, K., Aðalgeirsdóttir, G., and Cappelen, J.: Quantifying energy and mass fluxes controlling Godthåbsfjord freshwater input in a 5-km simulation (1991–2012), J. Climate, 28, 3694–3713, 2015. a

Langen, P. L., Fausto, R. S., Vandecrux, B., Mottram, R. H., and Box, J. E.: Liquid water flow and retention on the Greenland ice sheet in the regional climate model HIRHAM5: Local and large-scale impacts, Front. Earth Sci., 4, 110,, 2017. a

Lasher, G. E., Axford, Y., McFarlin, J. M., Kelly, M. A., Osterberg, E. C., and Berkelhammer, M. B.: Holocene temperatures and isotopes of precipitation in Northwest Greenland recorded in lacustrine organic materials, Quaternary Sci. Rev., 170, 45–55, 2017. a, b

Latif, M. and Keenlyside, N. S.: A perspective on decadal climate variability and predictability, Deep-Sea Res. Pt. II, 58, 1880–1894, 2011. a

Lecavalier, B. S., Milne, G. A., Simpson, M. J., Wake, L., Huybrechts, P., Tarasov, L., Kjeldsen, K. K., Funder, S., Long, A. J., Woodroffe, S., Dyke, A., and Larsen, N.: A model of Greenland ice sheet deglaciation constrained by observations of relative sea level and ice extent, Quaternary Sci. Rev., 102, 54–84, 2014. a

Lecavalier, B. S., Fisher, D. A., Milne, G. A., Vinther, B. M., Tarasov, L., Huybrechts, P., Lacelle, D., Main, B., Zheng, J., Bourgeois, J., and Dyke, A.: High Arctic Holocene temperature record from the Agassiz ice cap and Greenland ice sheet evolution, P. Natl. Acad. Sci. USA, 114, 5952–5957, 2017. a

Lesnek, A. J. and Briner, J. P.: Response of a land-terminating sector of the western Greenland Ice Sheet to early Holocene climate change: Evidence from 10Be dating in the Søndre Isortoq region, Quaternary Sci. Rev., 180, 145–156, 2018. a

Liu, Z., Otto-Bliesner, B., He, F., Brady, E., Tomas, R., Clark, P., Carlson, A., Lynch-Stieglitz, J., Curry, W., Brook, E., Erickson, D., Jacob, R., Kutzbach, J., and Cheng, J.: Transient simulation of last deglaciation with a new mechanism for Bølling-Allerød warming, Science, 325, 310–314, 2009. a, b, c, d, e

Liu, Z., Carlson, A. E., He, F., Brady, E. C., Otto-Bliesner, B. L., Briegleb, B. P., Wehrenberg, M., Clark, P. U., Wu, S., Cheng, J., Zhang, J., Noone, D., and Zhu, J.: Younger Dryas cooling and the Greenland climate response to CO2, P. Natl. Acad. Sci. USA, 109, 11101–11104, 2012. a

McFarlin, J. M., Axford, Y., Osburn, M. R., Kelly, M. A., Osterberg, E. C., and Farnsworth, L. B.: Pronounced summer warming in northwest Greenland during the Holocene and Last Interglacial, P. Natl. Acad. Sci. USA, 115, 6357–6362, 2018. a, b

Nash, J. and Sutcliffe, J.: River forcasting using conceptual models, 1. A discussion of principles, J. Hydrol., 10, 280–290, 1970. a

Nielsen, L. T., Aðalgeirsdóttir, G., Gkinis, V., Nuterman, R., and Hvidberg, C. S.: The effect of a Holocene climatic optimum on the evolution of the Greenland ice sheet during the last 10 kyr, J. Glaciol., 64, 477–488, 2018. a

Nusbaumer, J., Wong, T. E., Bardeen, C., and Noone, D.: Evaluating hydrological processes in the Community Atmosphere Model Version 5 (CAM5) using stable isotope ratios of water, J. Adv. Model. Earth Syst., 9, 949–977, 2017. a

Okazaki, A. and Yoshimura, K.: Development and evaluation of a system of proxy data assimilation for paleoclimate reconstruction, Clim. Past, 13, 379–393,, 2017. a

Okazaki, A. and Yoshimura, K.: Global evaluation of proxy system models for stable water isotopes with realistic atmospheric forcing, J. Geophys. Res.-Atmos., 124, 8972–8993, 2019. a

PAGES 2k-PMIP3 group: Continental-scale temperature variability in PMIP3 simulations and PAGES 2k regional temperature reconstructions over the past millennium, Clim. Past, 11, 1673–1699,, 2015. a

Pausata, F. S. and Löfverström, M.: On the enigmatic similarity in Greenland δ18O between the Oldest and Younger Dryas, Geophys. Res. Lett., 42, 10–470, 2015. a, b

Peltier, W.: Global glacial isostasy and the surface of the ice-age Earth: the ICE-5G (VM2) model and GRACE, Annu. Rev. Earth Pl. Sc., 32, 111–149, 2004. a

Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheet-shelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295,, 2012. a

Rasmussen, S. O., Andersen, K. K., Svensson, A., Steffensen, J. P., Vinther, B. M., Clausen, H. B., Siggaard-Andersen, M.-L., Johnsen, S. J., Larsen, L. B., Dahl-Jensen, D., Bigler, M., Röthlisberger, R., Fischer, H., Goto-Azuma, K., Hansson, M., and Ruth, U.: A new Greenland ice core chronology for the last glacial termination, J. Geophys. Res.-Atmos., 111, D06102,, 2006. a, b

Rasmussen, S. O., Abbott, P. M., Blunier, T., Bourne, A. J., Brook, E., Buchardt, S. L., Buizert, C., Chappellaz, J., Clausen, H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guillevic, M., Kipfstuhl, S., Laepple, T., Seierstad, I. K., Severinghaus, J. P., Steffensen, J. P., Stowasser, C., Svensson, A., Vallelonga, P., Vinther, B. M., Wilhelms, F., and Winstrup, M.: A first chronology for the North Greenland Eemian Ice Drilling (NEEM) ice core, Clim. Past, 9, 2713–2730,, 2013. a, b

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28, 2014. a

Raymond, C. F.: Deformation in the vicinity of ice divides, J. Glaciol., 29, 357–373, 1983. a

Reeh, N. and Fisher, D.: Noise in accumulation rate and δ18O time series as determined from comparison of adjacent Greenland and Devon Island ice cap cores, Ottawa, Polar Continental Shelf Project (Internal Report), 1983. a

Robin, G. d. Q.: Ice cores and climatic change, Philos. T. Roy. Soc. B, 280, 143–168, 1977. a

Roe, G. H. and Lindzen, R. S.: The mutual interaction between continental-scale ice sheets and atmospheric stationary waves, J. Climate, 14, 1450–1465, 2001. a

Roy, K. and Peltier, W.: Relative sea level in the Western Mediterranean basin: A regional test of the ICE-7G_NA (VM7) model and a constraint on Late Holocene Antarctic deglaciation, Quaternary Sci. Rev., 183, 76–87, 2018. a

Schüpbach, S., Fischer, H., Bigler, M., Erhardt, T., Gfeller, G., Leuenberger, D., Mini, O., Mulvaney, R., Abram, N. J., Fleet, L., Frey, M. M., Thomas, E., Svensson, A., Dahl-Jensen, D., Kettner, E., Kjaer, H., Seierstad, I., Steffensen, J. P., Rasmussen, S. O., Vallelonga, P., Winstrup, M., Wegner, A., Twarloh, B., Wolff, K., Schmidt, K., Goto-Azuma, K., Kuramoto, T., Hirabayashi, M., Uetake, J., Zheng, J., Bourgeois, J., Fisher, D., Zhiheng, D., Xiao, C., Legrand, M., Spolaor, A., Gabrieli, J., Barbante, C., Kang, J.-H., Hur, S. D., Hong, S. B., Hwang, H. J., Hong, S., Hansson, M., Iizuka, Y., Oyabu, I., Muscheler, R., Adolphi, F., Maselli, O., McConnell, J., and Wolff, E. W.: Greenland records of aerosol source and atmospheric lifetime changes from the Eemian to the Holocene, Nat. Commun., 9, 1476,, 2018. a

Severinghaus, J. P. and Brook, E. J.: Abrupt climate change at the end of the last glacial period inferred from trapped air in polar ice, Science, 286, 930–934, 1999. a, b, c

Severinghaus, J. P., Sowers, T., Brook, E. J., Alley, R. B., and Bender, M. L.: Timing of abrupt climate change at the end of the Younger Dryas interval from thermally fractionated gases in polar ice, Nature, 391, 141–146, 1998. a, b, c, d, e

Simpson, M. J., Milne, G. A., Huybrechts, P., and Long, A. J.: Calibrating a glaciological model of the Greenland ice sheet from the Last Glacial Maximum to present-day using field observations of relative sea level and ice extent, Quaternary Sci. Rev., 28, 1631–1657, 2009. a

Steig, E. J., Grootes, P. M., and Stuiver, M.: Seasonal precipitation timing and ice core records, Science, 266, 1885–1887, 1994. a, b

Steiger, N. J., Steig, E. J., Dee, S. G., Roe, G. H., and Hakim, G. J.: Climate reconstruction using data assimilation of water isotope ratios from ice cores, J. Geophys. Res.-Atmos., 122, 1545–1568, 2017. a

Stenni, B., Masson-Delmotte, V., Selmo, E., Oerter, H., Meyer, H., Röthlisberger, R., Jouzel, J., Cattani, O., Falourd, S., Fischer, H., Hoffmann, G., Iacumin, P., Johnsen, S., Minster, B., and Udisti, R.: The deuterium excess records of EPICA Dome C and Dronning Maud Land ice cores (East Antarctica), Quaternary Sci. Rev., 29, 146–159, 2010. a

Stuiver, M. and Grootes, P. M.: GISP2 oxygen isotope ratios, Quaternary Res., 53, 277–284, 2000. a, b

Svensson, A., Andersen, K. K., Bigler, M., Clausen, H. B., Dahl-Jensen, D., Davies, S. M., Johnsen, S. J., Muscheler, R., Rasmussen, S. O., Röthlisberger, R., Steffensen, J., and Vinther, B.: The Greenland ice core chronology 2005, 15–42 ka. Part 2: comparison to other records, Quaternary Sci. Rev., 25, 3258–3267, 2006. a, b

Teutschbein, C. and Seibert, J.: Bias correction of regional climate model simulations for hydrological climate-change impact studies: Review and evaluation of different methods, J. Hydrol., 456, 12–29, 2012. a

Vinther, B. M., Clausen, H. B., Johnsen, S. J., Rasmussen, S. O., Andersen, K. K., Buchardt, S. L., Dahl-Jensen, D., Seierstad, I. K., Siggaard-Andersen, M.-L., Steffensen, J. P., Svensson, A., Olsen, J., and Heinemeier, J.: A synchronized dating of three Greenland ice cores throughout the Holocene, J. Geophys. Res.-Atmos., 111, D13102,, 2006. a, b

Vinther, B. M., Clausen, H. B., Fisher, D., Koerner, R., Johnsen, S. J., Andersen, K. K., Dahl-Jensen, D., Rasmussen, S. O., Steffensen, J. P., and Svensson, A.: Synchronizing ice cores from the Renland and Agassiz ice caps to the Greenland Ice Core Chronology, J. Geophys. Res.-Atmos., 113, D08115,, 2008. a

Vinther, B. M., Buchardt, S. L., Clausen, H. B., Dahl-Jensen, D., Johnsen, S. J., Fisher, D., Koerner, R., Raynaud, D., Lipenkov, V., Andersen, K. K., Blunier, T., Rasmussen, S., Steffensen, J., and Svensson, A.: Holocene thinning of the Greenland ice sheet, Nature, 461, 385–388, 2009.  a, b, c, d, e

Vizcaino, M., Mikolajewicz, U., Ziemen, F., Rodehacke, C. B., Greve, R., and Van Den Broeke, M. R.: Coupled simulations of Greenland Ice Sheet and climate change up to AD 2300, Geophys. Res. Lett., 42, 3927–3935, 2015. a

Werner, M., Mikolajewicz, U., Heimann, M., and Hoffmann, G.: Borehole versus isotope temperatures on Greenland: Seasonality does matter, Geophys. Res. Lett., 27, 723–726, 2000. a, b

Whitaker, J. S. and Hamill, T. M.: Ensemble data assimilation without perturbed observations, Mon. Weather Rev., 130, 1913–1924, 2002. a

Young, N. E. and Briner, J. P.: Holocene evolution of the western Greenland Ice Sheet: Assessing geophysical ice-sheet models with geological reconstructions of ice-margin change, Quaternary Sci. Rev., 114, 1–17, 2015. a, b, c