Concurrent Changepoints in Greenland Ice Core δ¹⁸O Records and the North Atlantic Oscillation over the Past Millennium

Abstract

Structural changes, or changepoints, coinciding in multiple ice core records over the Greenland Ice Sheet (GrIS) may reflect a widespread response of the GrIS to atmospheric forcing. Thus, to better understand how atmospheric circulation may regulate sudden changes in δ¹⁸O of Greenland precipitation, we seek synchronous changepoints occurring in ice core-derived δ¹⁸O time series across the GrIS and in the North Atlantic Oscillation (NAO) over the past millennium. By utilizing a Bayesian changepoint detection method, four changepoint horizons were revealed: at the beginning of the 20th century, in the late-15th century, and around the turn of the 11th and 10th centuries. Although the changepoints in ice core δ¹⁸O records exhibited distinctive spatial arrangements in each horizon, all corresponded to changepoints in the NAO, indicative of a consistent atmospheric influence on GrIS surface changes over the past millennium.

1. Introduction

Atmospheric circulation overlying the Greenland Ice Sheet (GrIS) shows marked interannual-to-decadal oscillations arising from large-scale natural oceanic/atmospheric modes of variability [1,2,3] such as the North Atlantic Oscillation (NAO; [4]). The NAO is the leading mode of atmospheric circulation variability in the North Atlantic [4] described by an index (NAO index), most often calculated from sea level pressure, and is widely used as a general indicator for the strength of the westerlies over the eastern North Atlantic [5]. It is known to influence the location of storm tracks across the North Atlantic and to regulate Greenland surface conditions by mediating opposing atmospheric circulation patterns in its different phases [6,7,8]. NAO phase changes may occur gradually or relatively abruptly in time, which, either way, can be interpreted as a change in the distribution of the NAO-index values. However, the detection of such changepoints’ abundancy and their typical magnitudes is constrained by relatively short observational records.

Thus, placing observed coupled changes in GrIS surface conditions and atmospheric circulation aloft into the context of paleoclimatic time scales is necessary to better understand the consistency of their co-variability with regard to synchronous sudden changes. A vast number of high-resolution ice core-derived paleoclimate proxy parameters are available from the GrIS [9,10,11], some of which approximate ambient temperatures—e.g., firn/ice δ¹⁸O records [12]—that altogether offer insights into past Greenland changes beyond the observational era [13]. The water stable isotopes from firn and ice profiles are one of the most studied geochemical parameters when it comes to investigating past climates [12]. Firn/ice stable isotope signature is known to be affected by changes in evaporation conditions and air mass trajectories [14], frequency of weather regimes [15], and circulation patterns [16,17].

Studies of glaciochemical and water isotopic parameters on single ice cores have shown that if atmospheric circulation changes, then the change in moisture source region [18] and transport path of precipitation [19] is rapidly reflected in the water stable isotopic composition of fallen snow. A small set of sea salt aerosol proxy records from the northern GrIS show a consistent relationship with large-scale pressure patterns and storm activity over the North Atlantic region [10]. However, the characteristics of this pattern could have not have been reconstructed from a single ice core due to the low signal-to-noise ratio. A dominant influence of the NAO was found on the west-central GrIS in the early 20th century by studying an extensive collection of ice core-derived annual accumulation histories [20]. The leading mode of variability of ice core δ¹⁸O for 1778–1970 CE of multiple ice cores (n = 14) was linked to the NAO primarily in the Southern region of the GrIS [11]. A subsequent study, spanning a slightly broader time interval compared to [11] and reinforced by isotope simulations, evaluated an extended set (n = 22) of shallow ice core derived δ¹⁸O series across the GrIS and pointed out spatially variable depletion/enrichment in the isotopic composition of ice/firn depending on the dominant NAO mode [15].

In an effort to complement previous studies, ice core-derived water stable isotope records available from the GrIS and its vicinity are explored to find concurrent changepoints that not only occur in multiple ice cores but are also coincident with structural changes in the NAO. Major centennial climate anomalies encompassed by the past millennium, such as the Medieval Climate Anomaly (from ~950 to ~1250 CE) or the Little Ice Age (from ~1450 to ~1850 CE) [21], might invite to query whether synchronous changes in ice core δ¹⁸O and the NAO reflect those periods in the past millennium.

If multiple objectively defined and statistically significant changes occur among different ice cores, it would clearly suggest that these oscillations are not merely artifacts of local variability of the GrIS but instead are triggered by larger-scale forcing(s) [22]. Hence, the present study complements previous efforts and elucidates the persistence of the observed strong physical link between precipitation δ¹⁸O variability and atmospheric circulation beyond the instrumental era. Although the underlying mechanisms causing abrupt Greenland climate changes are far from fully understood, the goal of this study is to explore statistical relationships instead of seeking driving mechanisms quantitatively.

2. Materials and Methods

Ice core-derived δ¹⁸O records were obtained from the Iso2k database [23,24] utilizing all ice core-derived δ¹⁸O records with annual resolution poleward of 59° N from the GrIS and its vicinity (n = 30; Figure 1 and Figure S3). This improved spatial coverage compared with previous multi-ice core studies—e.g., Vinther et al. [11] and Ortega et al. [15]—helps to identify possible regional imprints inherent to the synoptic patterns influencing Greenland climate variability.

Ice core δ¹⁸O record availability steadily decreases from the 19th century backward (Figure 1). The five records available before 750 CE represent only the islands of the Canadian Arctic and the western part of Greenland. Therefore, our investigation is restricted to the past millennium, where available records ($n\ge 10$) provided complete coverage for the study area (Figure 1) that overlap with instrumental climate observations with the density of ice core records being maximal, permitting to take advantage of the spatio-temporal structure of the isotopic variability [15]. The extended instrumental winter (DJF) North Atlantic Oscillation record (1865–2002 CE; [25]) calculated as the difference in mean monthly surface pressure anomalies between Iceland and the Azores was used as an observational benchmark hereafter referred to as NAO.

Regarding the pre-instrumental era, there are multiple NAO reconstructions available based on different methodologies and predictors; however, reconstructions relying on single or “too few” proxy locations have a high risk of insufficiently describing the NAO variability [26]. Therefore, in the present study, only those reconstructions were considered, which used an extensive proxy network and fully covered the past millennium, sufficiently represented by the studied ice core dataset. This criterion yielded two NAO reconstructions that represent large-scale modes of atmospheric variability aloft the ice sheet: a model-tested NAO reconstruction [27] covering 1049–1969 CE (NAO*), and a winter NAO index using a network of 97 Euro–Mediterranean tree-ring series covering 910–2018 CE (NAO_med; [28]). The NAO is significantly correlated with the Greenland blocking index [11,29] calculated from 500 hPa geopotential heights. Thus, for simplicity, the analyses are restricted to the NAO and consider it as a proxy for changes in atmospheric forcing on the ice sheet while noting that similar analysis may be performed using other climate variability indices.

Since some of the records lacked high-frequency variability in the sub-decadal range (Figure S3), a 10-year lowpass filter was applied to all the ice core records to achieve common variability. If this had not been conducted, the different spectral characteristics of the original ice core δ¹⁸O records would have required different parameter settings causing unnecessary difficulties. The filtering also increases the signal-to-noise ratio relative to the magnitude of the supposed changes, which aids the detectability of valid changepoints [30,31]. The same lowpass filter was applied to the other indices for consistency. The 10-yr low-passed time series were then centered and subjected to the linear model of the Bayesian changepoint algorithm of Ruggieri [32] (BCPa for short), an approach that can generate the posterior distribution on both the number and location of changepoints in a data set.

The changepoint analysis has been an active area of statistical research for many years, e.g., Scott and Knott [33]. In every case, the goal is to perform a statistical test to determine if a homogeneous model or a piecewise model provides a more appropriate fit to the data, and if piecewise, how many breaks are appropriate. If a piecewise model provides the more appropriate fit, then the point at which the model changes is known as a “changepoint”. The actual model used in each approach varies by application ranging from a constant model (which can be used to detect changes in the mean signal) to a more general regression model.

Applications of changepoint detection in a climate setting are numerous and include NAO [34], carbon dioxide concentrations [35], tropical cyclones [36], precipitation [37], storm intensity [38], temperature proxy records ([39]: Sect. 4.1. therein) and global surface temperatures [32], among others. The literature on changepoint models (generally) and their applications to climate data are vast, so a full review is beyond the scope of this paper.

We define the point at which the statistical properties of a model change as a “changepoint”. For example, when using a linear (e.g., trend) model, $Y={\beta }_0+{\beta }_{1}t+ϵ,$ any change to the slope (${\beta }_1)$, intercept (${\beta }_0$), or variance (${\sigma }^2$) in the time series indicates the presence of a changepoint. Ultimately, the goal is to fit a piecewise regression model to each time series. BCPa assumes that the parameters of the model for any two segments of the data are independent (i.e., a product partition model [40]). In this context, each segment of the data represents a different climate regime, and a “changepoint” represents the boundary between.

BCPa has three steps:

• Calculate the probability of the data given the model between any two time points, P(i,j). This represents all possible climate regimes (in terms of timing) in the data set;
• Using the probabilistic calculations from step (1) as building blocks, recursively piece together these segments, adding one changepoint at a time. Specifically, if we define P_k(j) to be the probability that the first j data points contain k changepoints, then

  $$P_k\left(j\right)={\displaystyle \sum }_{v<j}{P}_{k-1}\left(v\right)P\left(v+1,j\right)$$

  for k = 1 to k_max, where P(v + 1, j) was calculated in step (1);
• Using Bayes’ rule, sample from the posterior distribution on the number and location of the changepoints, as well as the parameters of the model in each segment.

The specific parameters applied were as follows: minimum distance between the possible changepoints d_min= 10 steps (years in the present case), hyper parameter for the prior on the regression coefficients k₀ = (0.001, 0.01), prior on the variance of the residuals, v₀ = 1, σ₀ = variance of each time series, the maximum number of changepoints k_max = 50, and the number of sampled solutions num.samp = 10,000.

The reader interested in a general overview of changepoint analysis is referred to Ruggieri [41], while those interested in a more technical derivation of BCPa are referred to Ruggieri [32] and SOM Sect. 1 (Figure S1).Those interested in the description of the specific parameters are referred to SOM Sect. 2 (Figure S2), while the software used is described in SOM Sect. 3.

A Bayesian approach to the changepoint problem provides a posterior distribution, i.e., uncertainty estimates, on the location of changepoints. Therefore, the obtained changepoint time series were passed through a 5-yr centralized rolling summary function in order to cumulate the probability dispersed between consecutive years. The use of this function assumes that the uncertainty in the location of a changepoint is less than 5 years. Allowing for a larger window could increase the overall probability of a changepoint (depending on how abrupt and how “obvious” the change is) in a certain region of the data. This cumulative probability is the probability of the changepoint(s) shown in the study. Only “practically certain” changepoints with a 5-yr sum probability >99% were considered in the evaluation.

In addition, since “paleoData_chronologyIntegrationTimeUncertainty” was determined as between “1 and 10” years in > 90% of the number of assessed ice cores, a 10-year window was considered when searching for coinciding changepoints. These coinciding changepoints were taken as a “changepoint horizon” if at least ~25% of the ice cores available at that time indicated a changepoint. Here, the assumption is that if changepoints are correlated, then they should occur at a (relatively) similar point in time. For example, if a changepoint centered in one ice core in the year 1921 with a 99.4% probability and a different ice core in 1929 with a 99.1%, then these were considered to fall into the same so-called changepoint horizon. It was also investigated how the ice cores representative of a changepoint horizon are clustered in space.

3. Results and Discussion

The main goal of the study was to reveal how concurrent changes in Greenland surface conditions and the NAO may have occurred synchronously over and beyond the instrumental era.

We began by describing the timing of the events associated with synchronous changes in the GrIS ice core δ¹⁸O records. These constitute four changepoint horizons over the last millennium centered at years 1933 (Figure 2), 1460, 1190, and 1085 CE (Figure 3) referred to hereinafter as CpH₁₉₃₃, CpH₁₄₆₀, CpH₁₁₉₀, and CpH₁₀₈₅, respectively. Out of these, the relatively most abundant was the CpH₁₁₉₀ with 40% of the available ice cores (n_total = 10), indicating a changepoint at > 99% probability. Another 20% indicated a changepoint at >85% probability. Meanwhile, CpH₁₉₃₃ was indicated by the greatest number of ice cores (seven out of n_total = 30). Five additional ice core δ¹⁸O records show a changepoint between 1928 and 1937 CE with a posterior probability > 90%, reinforcing the presence and spatial representativity of the CpH₁₉₃₃. These results are robust in terms of the chosen lowpass filter, as the same set of ice core δ¹⁸O changepoints were found using 13-yr filtering; 78% with the same probability (>99%) as in the case of the 10-yr lowpass, and 22% with a slightly lower probability (~95%).

The individual ice cores defining a changepoint horizon were clustered in space. These reflect the three major regions of the GrIS delineated based on a synoptic survey of precipitation and possible effects of orography on moving cyclones [42]. This regional division is used to describe the spatial pattern(s) of the CpHs. The most recent changepoint horizon found in the second quarter of the 20th century (CpH₁₉₃₃) is well represented by more than half of the ice core δ¹⁸O records from the Southern Region (Figure 2). The horizon from the mid-15th century (CpH₁₄₆₀) covers the North Central Region together with an additional ice core δ¹⁸O record from the Northern region (B21) (Figure 3b). Meanwhile, in the oldest horizons, CpH₁₁₉₀ and CpH₁₀₈₅, ice cores from the Central region were involved, accompanied by a core from the Ellesmere Island (A79; Figure 3c,d).

A single ice core (B19) located in the NE sector of the GrIS can be considered the “most responsive”, since it is in all CpHs (Figure 2b and Figure 3b–d). This is in line with the expectation that precipitation δ¹⁸O in this area is relatively more sensitive to changes in moisture uptake conditions compared with other parts of the island [43].

3.1. Common Changepoint in Ice Core δ¹⁸O and Atmospheric Circulation in the 1930s

The most recent changepoint horizon is centered in the year 1933 CE (CpH₁₉₃₃, Figure 2) and is the only one found since the beginning of the industrial era. This changepoint corresponds to the well-known sharp regime shift in atmospheric circulation across the Atlantic–Arctic sector, indicating the end of a generally colder era in Greenland [44]. This, in the more recent literature, is referred to as the “early 20th century warming” [45,46,47]. Taken together, concurrent changes in both the NAO and NAO* and in the ice core δ¹⁸O records mark a transition towards generally warmer Greenland conditions starting from the early 20th century up until the 1960s (Figure 2a). As for the ice cores, the 1933 horizon is defined by seven (probability > 99%) and five (probability > 90%) out of 30 cores across the GrIS (Figure 2b).

At CpH₁₉₃₃, the ice core δ¹⁸O records from the Southern Region show a unanimous positive shift, while the NAO shifts greatly from positive to negative (Figure 2). This is in harmony with simulations derived from an isotope-enabled general circulation model, where the negative NAO phase is accompanied by positive δ¹⁸O anomalies characterizing the complete southern region and West Central Greenland [15], and can be explained by the NAO’s influence on storm tracks [4,48]; thus on the moisture source regions and distillation pathways. Such fluctuations are expected to be reflected in the firn/ice δ¹⁸O signal representing the different sectors of the GrIS [43,49].

The dominant moisture source location is known to vary with the NAO phase across the GrIS most strongly for the northeastern and southeastern sectors of the GrIS. During the NAO⁺ phase, uptake locations are shifted to the Arctic region (north of 60° N), while in NAO^– phase, more long-range moisture transport from the Atlantic (south of 50° N) was identified [43]. However, different regions of the GrIS experience distinct changes in the average moisture source locations, which is manifested in the stable isotopic signal of ice/firn. The strength and direction of this relationship vary across Greenland and its surroundings [15].

3.2. Common Changepoints in Ice Core δ¹⁸O and the NAO over the Past Millennium

Changepoint horizons prior to the instrumental era are also explored to elucidate the persistency of coupled changes in Greenland ice cores and the NAO. The mid-15th century horizon (CpH₁₄₆₀) was formed by 25% of the available ice core δ¹⁸O records (n = 4) (Figure 3b). At CpH₁₄₆₀, NAO* shows a recovery from a decadal negative excursion (Figure 3a). In the North and North Central Region of the GrIS, where the ice cores indicating the CpH₁₄₆₀ are located (Figure 3b), the moisture source region shifts from more distant to closer as a result of the NAO changing phase from negative to positive [43]. The timing of the CpH₁₄₆₀ matches the decadal shift in broad-scale moisture variability on both sides of the Atlantic, a shift that has also been linked to a negative NAO anomaly in the middle of the 15th century [50] associated with the start of the Little Ice Age. Moreover, this corresponds to a regime shift of atmospheric processes at centennial time scales, i.e., to the transition from a more probable NAO^– epoch (1500–1800 CE) to a more probable NAO⁺ epoch (1200–1400 CE) [51]. Furthermore, CpH₁₄₆₀ is spatially distinct from CpH₁₉₃₃, as it is the only one showing changepoints in the northern GrIS.

In the penultimate documented horizon (CpH₁₁₉₀), 40% of the available ice core δ¹⁸O records show a changepoint at a rather common negative peak, which is followed by a multidecadal period of less variability with a relatively higher mean (Figure 3c). In the case of the earliest horizon (CpH₁₀₈₅), the time series of ice core δ¹⁸O show a common negative shift (Figure 3d). The two CpHs correspond closely to the modeled change in the multidecadal probabilistic distribution of NAO modes. The two pre-13th century changepoint horizons (CpH_{1190 and 1085}) were also mirrored in the NAO_med reconstruction. The timing of the CpH₁₁₉₀ is in harmony with the change from more probable NAO⁺ epoch (1200–1400 CE) to more probable NAO^– epoch (1100–1200 CE), and CpH₁₀₈₅ matches the transition from a more probable NAO⁻ epoch, i.e., 1100–1200 CE, to a more probable NAO⁺ epoch, i.e., 1000–1100 CE [51].

The CpH₁₁₉₀ (Figure 3a,c) shows an opposite NAO phase change compared with the CpH₁₄₆₀ (Figure 3a,b), and hence the corresponding coupled changes in atmospheric circulation imprinted in ice core water stable isotopes [15,17,52] are likely to feature similarly contrasting patterns between the two CpHs.

The characteristic difference between the industrial and pre-industrial era changepoint is manifested in the spatial clustering of the corresponding ice cores in, basically, the south and north part of the GrIS. This divide is also mirrored in the dichotomic character of the modeled stable isotope maps of Greenland [15].

It should be noted that there were three additional changepoints detected in the NAO* time series besides the previously mentioned ones. One of these was found in the mid-19th century and two in the 12th century. The oldest of these, occurring in the mid-12th century, is the most noteworthy since it coincides with changepoints in the NAO (1147 CE) and both the Crete and GISP2 records (posterior probability > 99.5%), which are in relatively close vicinity ~150 km apart (Figure 1a). Although only two out of the ten available ice cores indicated a changepoint, which does not meet the 25% criteria of the study to be considered a CpH, this finding reinforces the systematic and spatially clustered response of Greenland ice core δ¹⁸O records to changes in the NAO. The proportion of ice cores showing a changepoint was lower for the other two NAO* changepoints, but still greater than zero in both cases.

4. Summary

This study shows the first joint assessment for detecting changepoint horizons in the ice core derived water stable isotope records with annual resolution available from the GrIS and its vicinity for the past millennium, and also explores if any coinciding abrupt changes exist in the dominant mode of large-scale atmospheric circulation in the North Atlantic Basin. The four detected changepoint horizons (centered at years 1933, 1460, 1190, and 1085 CE) were accompanied by changepoints in the North Atlantic Oscillation, confirming the hypothesized consistency of the link between large-scale atmospheric circulation regimes and Greenland ice/firn δ¹⁸O on a decadal time scale. The 20th century changepoint horizon detected in the Greenland ice/firn δ¹⁸O records was reflected in both the instrumental and the reconstructed NAO index time series, while the pre-instrumental changepoint horizons corresponded to the explicit change in the multidecadal probabilistic distribution of NAO [51] modes lending further credibility to the presented results.

Moreover, the spatial clustering of the ice core δ¹⁸O records indicates that a given changepoint horizon was found to have an encouraging consensus with subregions of the GrIS delineated with (isotope enabled) moisture source models. Specifically, the E–NE areas, which are recognized in the literature as the most sensitive and responsive to changes in moisture source regions, were found to be represented in all changepoint horizons.