Impact of abrupt sea ice loss on Greenland water isotopes during the last glacial period

Significance The Dansgaard–Oeschger events contained in Greenland ice cores constitute the archetypal record of abrupt climate change. An accurate understanding of these events hinges on interpretation of Greenland records of oxygen and nitrogen isotopes. We present here the important results from a suite of modeled Dansgaard–Oeschger events. These simulations show that the change in oxygen isotope per degree of warming becomes smaller during larger events. Abrupt reductions in sea ice also emerge as a strong control on ice core oxygen isotopes because of the influence on both the moisture source and the regional temperature increase. This work confirms the significance of sea ice for past abrupt warming events.

Greenland ice cores provide excellent evidence of past abrupt climate changes. However, there is no universally accepted theory of how and why these Dansgaard-Oeschger (DO) events occur. Several mechanisms have been proposed to explain DO events, including sea ice, ice shelf buildup, ice sheets, atmospheric circulation, and meltwater changes. DO event temperature reconstructions depend on the stable water isotope (δ 18 O) and nitrogen isotope measurements from Greenland ice cores: interpretation of these measurements holds the key to understanding the nature of DO events. Here, we demonstrate the primary importance of sea ice as a control on Greenland ice core δ 18 O: 95% of the variability in δ 18 O in southern Greenland is explained by DO event sea ice changes. Our suite of DO events, simulated using a general circulation model, accurately captures the amplitude of δ 18 O enrichment during the abrupt DO event onsets. Simulated geographical variability is broadly consistent with available ice core evidence. We find an hitherto unknown sensitivity of the δ 18 O paleothermometer to the magnitude of DO event temperature increase: the change in δ 18 O per Kelvin temperature increase reduces with DO event amplitude. We show that this effect is controlled by precipitation seasonality.
There has recently been significant progress in reconstructing abrupt DO temperatures increases over Greenland using nitrogen isotopes δ 15 − N2 (12,19). This work indicates jumps in temperature over Greenland of up to 16.5 ± 3 K within a few decades (12,19). A logical but challenging next step is to elucidate how geographical patterns of change in key ice core records, particularly δ 18 O, from Greenland ice cores can be used to provide that crucial missing information on the nature and cause of abrupt warming events, sea ice loss, and its relationship to abrupt temperature rises (20,21).
DO events are imprinted across the whole of Greenland: wherever last glacial ice is preserved, ice core measurements capture these events (10-12, 19, 22). However, the magnitude of the DO imprint is not identical across the Greenland ice sheet. Early DYE3 ice core measurements suggest that δ 18 O changes during DO warmings may be larger in the south of Greenland (10) compared with central Greenland. More recent ice core data ( Fig. 1  A and B) imply that, while the magnitudes of DO temperature and accumulation changes (from δ 15 − N2 and δ 18 O) are larger in central Greenland compared with the north and northwest (12), δ 18 O changes could be larger and are likely more variable in the north and northwest compared with central sites (9,12,22). How this spatial variability relates to sea ice loss is currently unknown.
General circulation model (GCM) simulations of δ 18 O enable robust interpretation of records recovered from Greenland ice cores. In particular, they allow us to probe influences on the geographical patterns on the measured δ 18 O change. The ability to decode DO changes from δ 18 O records from Greenland ice cores is thus vital to test ideas about drivers of past abrupt climate change (20,(23)(24)(25). Here, we present results from a large ensemble of 32 isotope-enabled GCM simulations of DO-type events.
Our DO-type simulations use a freshwater hosing-type setup. Salt is progressively lost from the North Atlantic during stadial periods; classic hosing (26,27) mechanisms explain the stadial North Atlantic region cooling that we simulate. After a switch of forcing in the hosing, salt returns to the North Atlantic from the tropical Atlantic Ocean and the wider global ocean. This causes the onset of an abrupt warming DO event.
We generate a suite of stadial climates from a 1,500-y glacial period spin-up simulation and then branch 32 different simulations of DO-type warming events off from this range of stadial climates ( Fig. 1 C-E and SI Appendix, Fig. S1). The simulations feature different amplitudes of effective salt fluxes alongside the range of initial stadial sea ice states. When calculating stadialinterstadial differences, 50 y of data are used for each climate: the 20 y on either side of the salt flux switch are excluded.

Significance
The Dansgaard-Oeschger events contained in Greenland ice cores constitute the archetypal record of abrupt climate change. An accurate understanding of these events hinges on interpretation of Greenland records of oxygen and nitrogen isotopes. We present here the important results from a suite of modeled Dansgaard-Oeschger events. These simulations show that the change in oxygen isotope per degree of warming becomes smaller during larger events. Abrupt reductions in sea ice also emerge as a strong control on ice core oxygen isotopes because of the influence on both the moisture source and the regional temperature increase. This work confirms the significance of sea ice for past abrupt warming events.   These simulations provide a means to unlock the Greenland δ 18 O records of abrupt DO climate change.
We compare our simulations with high-resolution isotopic records from the last glacial period available from the NGRIP, GRIP, and GISP2 ice cores (11,28,29). In addition to the available Greenland ice core δ 18 O data, a recent identification of DO events is used to locate the abrupt warming transitions for each event (30).

Comparison with Records of δ 18 O from Ice Cores
Comparison of simulated DO-type warmings in simulations with significant δ 18 O jumps with equivalent Greenland ice core measurements (of DO-type abrupt temperature rises) shows good model-data agreement in the magnitude and rate of the abrupt rises in δ 18 O ( Fig. 1 and SI Appendix, Fig. S1). Within uncertainties, the magnitudes of modeled precipitation and temperature increases (Table 1) are also in agreement with ice core observations (12,19). Increases in temperature and δ 18 O during DO events are always largest in the far south, with smaller changes in the north (Figs. 2A and 3). However, the modeled elevations are too low in the south (Materials and Methods). This likely contributes to the larger δ 18 O response at DYE3. A small region in northeast Greenland shows negative δ 18 O changes in some but not all simulations (Fig. 3).
Much of the near-Greenland sea ice loss tends to occur in the southwest on the east side of the Labrador Sea ( Fig. 2A). The geographical variability in δ 18 O individual simulations corresponds to the individual pattern of near-Greenland sea ice  loss ( Fig. 3 A-C). A similar geographical pattern also occurs in the temperature change patterns recorded in ice cores (12,20),and in modeled stadial-interstadial temperature changes (15), lending credence to these simulations. The δ 18 O increase is always strongly positive in the southwest of Greenland, including at DYE3, most commonly with a gradual reduction toward zero change in the northeast. However, some simulated DO events exhibit a stronger east-west geographical gradient in δ 18 O change (Fig. 3B), and others exhibit a stronger south-north gradient (Fig. 3C). Thus, although all simulations have large δ 18 O increases at the most southern DYE3 ice core site, δ 18 O increases are more variable between simulations at the central and northern sites ( Fig. 3 A-C). During some simulated events, NEEM changes are larger than those at GISP2 and GRIP and of the same magnitude of those at NGRIP.

Greenland Ice Core Paleothermometers
The δ 18 O-temperature, or traditional paleothermometer (31), coefficient (δ 18 O per Kelvin) over the DO warmings varies considerably across Greenland ( Fig. 2 B-D). To reconstruct DO temperature changes from δ 15 − N2, a model of firn densification and heat diffusion into the ice is required. This uses initial estimates of accumulation rate, ice age, and temperature, where the latter is derived from δ 18 O. Better independent constraints on accumulation rate and temperature from δ 18 O can thus reduce uncertainties on the reconstruction of abrupt DO temperature and accumulation change from δ 15 − N2 (12,19). It has previously been proposed that the δ 18 O-temperature relationship is dependent on obliquity insolation forcing (19). Here, we investigate another possibility. We calculate the paleothermometer coefficient associated with each individual simulated DO event (i.e., using coefficients from ∆δ 18 O and ∆ temperature) at each ice core site. With this approach applied to 15 simulations with significant (+2.0‰) DO rises in δ 18 O (SI Appendix, Fig. S1), we obtain a mean paleothermometer coefficient of 0.63 at DYE3; GISP2 yields 0.31. At GRIP, it is 0.30. NGRIP is 0.29, and at NEEM, it is 0.23‰ per Kelvin. A very similar geographical pattern emerges regardless the type of approach used to calculate the paleothermometer, with high values in the south and low values in the north. Within uncer-tainties, the simulated coefficients match the few observed coefficients (Materials and Methods) (12,19). Results here support the idea of considerable variability between DO events in the δ 18 O-temperature relationship (12,19).
Our simulations allow us to decipher what controls δ 18 Otemperature variability between DO events and between ice core sites. In our ensemble, the magnitude of the DO temperature increase has a strong control over the paleothermometer coefficient. Lower δ 18 O-temperature coefficients occur during larger DO events (Fig. 2B) (i.e., we find a strong systematic relationship between the size of the abrupt warming and the paleothermometer coefficient at all Greenland ice core sites). This finding also provides support for the idea that the paleothermometer is fundamentally dependent on the change in temperature at high latitude (32). The pattern varies over Greenland: small decreases occur in the paleothermometer coefficient with warming event size at DYE3 compared with large decreases in the coefficient at NEEM. This identification of a systematic dependence of the δ 18 O-temperature relationship on the size of the abrupt warming may be useful in further constraining Greenland abrupt temperature change records based on δ 15 − N2 and δ 18 O measurements.

Precipitation and δ 18 O Seasonality Changes
The current prevailing hypothesis is that most Greenland geographical variability in ∆δ 18 O is due to geographical variability in the seasonality of precipitation (9,12,19). To test this and to better understand the sea ice imprint of DO and the δ 18 O-temperature relationship, we calculate the impact of δ 18 O by only archiving climate information during periods of snow accumulation (33).
If a higher fraction of the yearly precipitation falls during winter months, the δ 18 O record will have a negative δ 18 O bias. In addition to this, the average δ 18 O in each month will also change as the climate moves from a stadial to an interstadial state. We quantify these two effects (by isolating the impact of changes in  The seasonal cycle of precipitation and δ 18 O both change during a DO event; a larger proportion of precipitation falls during colder months under the warmer interstadial climate relative to the cooler stadial climate. While these changes are less important in driving the majority of the geographical variability (or intercore differences) in δ 18 O across Greenland, compared with the pattern of near-Greenland sea ice loss, they are critical for understanding why sea ice controls on ∆δ 18 O vary so strongly across Greenland. In particular, the huge decreases in paleothermometer coefficients during larger warming events are dependent on changes in seasonality. Intensifications of precipitation seasonality under larger DO events reduce ∆δ 18 O everywhere in Greenland, but due to smaller values of ∆δ 18 Oseas in the north, the impact of precipitation seasonality is key in this region (Figs. 2 C and D and 3).
Intensifications of wintertime precipitation due to a larger wintertime area of open water around Greenland occur between stadials and interstadials (SI Appendix, Fig. S5). This registers as negative ∆Pseas (Fig. 3E and Table 1). Average ice core ∆Pseas is between −4.3 and −5.8‰, which exceeds the size of the recorded DO δ 18 O rise for four of five ice core sites.
Countering this, ∆δ 18 Oseas is considerably enriched, with an increase of between +4.2 and +12‰ across the ice core sites (Fig. 3F and Table 1). The change in evaporation in the ensemble is linearly dependent on sea ice coverage (SI Appendix, Fig. S6), with a strong dependence on the location where sea ice is reduced (SI Appendix, Fig. S7). This increase in (local) evaporation provides explanation for why accumulation and δ 18 O tend to rise in the vicinity of sea ice loss. As sea ice retreats during the interstadial, evaporation occurs much closer to Greenland (SI Appendix, Fig. S7); this moisture travels a much shorter distance to Greenland and therefore, is much less depleted than the moisture arriving over Greenland during the stadial when the sea ice edge is up to 20 • farther south.
∆Pseas is rather uniform across Greenland, but the geographical pattern of ∆δ 18 Oseas exerts considerable influence on how sea ice and temperature changes are recorded across Greenland. Thus, contrary to the reasonable current hypothesis that most Greenland geographical variability in ∆δ 18 O is due to geographical variability in the seasonality of precipitation (9,12,19), instead it is ∆δ 18 Oseas, which is much more important than ∆Pseas in driving the majority of geographical variability in ∆δ 18 O across Greenland.
While the huge increase in ∆δ 18 Oseas at the ice core sites more than compensates for negative ∆Pseas in the south, there is a fine balance in the central and northern regions between ∆Pseas and ∆δ 18 Oseas. This is due to the weak ∆δ 18 Oseas imprint. This results in weak relationships between δ 18 O and sea ice and the associated δ 18 O-temperature paleothermometer relationship (Fig. 2  B and C). Understanding these two large opposing changes, with distinctly different geographical patterns, is key to understanding Greenland DO δ 18 O changes.

Role of Sea Ice
In the north of Greenland, there is an extraordinary simulated range of DO event behaviors. This can be seen in the NEEM δ 18 O changes and the δ 18 O-temperature relationships at NEEM and partly, NGRIP ( Fig. 2 and Table 1): NEEM shows the largest range of δ 18 O-temperature coefficients (Fig. 2C), compared with any of the other ice cores sites. The amount of variance in δ 18 O for the largest simulated DO events directly explained by temperature (sea ice changes) is less than 29% (39%) at NEEM. However, at the other end of Greenland, the record of δ 18 O change at DYE3 for the largest events is nearly entirely 95% explained by sea ice changes (Fig. 2 C and D, r values converted to explained variances).
Although it is not possible to unambiguously attribute δ 18 O changes to particular components, like sea ice, temperature, atmospheric circulation, or storm tracks, the similar patterns of δ 18 O-temperature and δ 18 O-sea ice relationships (Fig. 2 C  and D), higher explained variances for sea ice over temperature (95% for sea ice vs. 92% for temperature at DYE3 and 70 vs. 62% at NGRIP), and sea ice impacts on the moisture sources and transport to Greenland suggest that sea ice exerts an even greater control on the stadial-interstadial δ 18 O over Greenland than temperature. This is because sea ice change controls both temperature in the wider region and the moisture availability.
This demonstration of the importance of the sea ice imprint on DO event Greenland ∆δ 18 O should help open the door to quantitative reconstructions of abrupt DO sea ice change based on these ice core measurements. Finally, these results show that precise and well-dated measurements of δ 18 O from DYE3 or from other new southern dome cores alongside careful isotopicenabled modeling (35,36) would be invaluable in quantifying changes in Arctic Sea ice during DO events.

Materials and Methods
Model Simulations. Global modeling studies with intermediate complexity models have shown bistability of the Atlantic Meridional Overturning Circulation (AMOC) in response to oscillatory or stochastic freshwater forcing in the North Atlantic (26). However, this behavior does not always appear in fully coupled 3D GCMs (18,37,38). Recently, however, one such model (Community Earth System Model 1) has been shown to predict significant nonforced oscillations in the strength of the AMOC and North Atlantic temperatures when subject to last glacial maximum (LGM) boundary conditions (27). This behavior has been termed a "kicked" salt oscillator. At the start of a DO-type warming, a large northward North Atlantic flux of salt occurs. This is associated with abruptly increased northward delivery of ocean heat and the melt of substantial areas of North Atlantic sea ice (5,6). A gradual external forcing related to orographic change of the North American Laurentide ice sheet or a small freshwater perturbation can also trigger similar oscillations (4,15). A version of the GCM called GFDL CM2Mc also exhibits nonforced oscillations of the AMOC, although only for a certain combination of background climate conditions (39). Here, we use a forced salt (or hosing) oscillation approach for setting up our suite of DO-type simulations.
Simulations are set up using an isotope-enabled version of the Hadley Center Coupled Model general circulation model (HadCM3). This GCM consists of a coupled atmosphere, ocean, and sea ice model and has been widely used to study past, present, and future climates (40,41). HadCM3 can also be run for hundreds of years on modern supercomputers (24,42,43). The ocean component of HadCM3 is a rigid lid model, with a fixed volume and water conservation through a time-invariant surface salinity flux that represents iceberg calving. Ref. 44 has the implementation of water isotope code in HadCM3. Ice sheets and sea ice in the model are initialized with δ 18 O values of −40 and −2‰, respectively. Using this model, salinity fluxes are applied in opposing directions to the North Atlantic vs. the rest of the global ocean surface. Simulations are set up using LGM ice sheets, orbital forcing, and greenhouse gas composition; additional details are in ref. 42. Every DO simulation is continued from the same standard 1,000-y glacial period (LGM) spin-up simulation. In the first instance, each DO-type simulation is then run for 500 y using an identical stadialtype forcing: a negative salinity flux, equivalent to 0.25 Sv of freshwater, is applied across the North Atlantic between 50  Fig. S1). Later interstadial negative salinity fluxes are applied using the same range of magnitudes. Each of these stadial-type simulations is branched off from a different year of the spin-up simulation. For each of the salt flux experiments, the surface layer of the ocean is freshened across the North Atlantic at a rate equivalent to the addition of between 0.1 and 1 Sv of freshwater, while the rest of the surface layer of the ocean is salinified at a rate equivalent to the loss of between 0.1 and 1 Sv (i.e., from −0.1 to −1 Sv) of freshwater. The suite of simulations is run using this stadialtype forcing for between 100 and 500 y, yielding a wide range of stadial climates. In each case, the North Atlantic (negative) salt forcing is always balanced by an equivalent-sized (positive) salt forcing applied across the rest of the global ocean; therefore, the net global freshwater (or salt) flux is always zero.
To generate a switch between stadial-type and interstadial-type climate states, a reversed salinity forcing, with positive salinity forcing over the North Atlantic and negative values over the remaining global ocean, is then applied. All forcings are between ±0.1 and ±1 Sv. The net global salinity flux is always zero. When we calculate stadial-interstadial difference, we use 50 y of data in each case: the 20 y on either side of the salt flux switch are excluded.
Greenland is represented by 80 grid points in HadCM3; thus, some smoothing of the surface topography is required to run the simulation (43). The modeled surface elevation is thus lower than the observed elevation. For the northern ice core sites, our modeled surface elevations are generally within 500 m of the present day surface elevations: NEEM: 2,450 vs. 2,341 m above sea level (observed vs. modeled, respectively). Similarly, NGRIP is 2,917 vs. 2,788 m above sea level (observed vs. modeled, respectively), and GISP2 is 3,216 vs. 2,832 m above sea level (observed vs. modeled, respectively). Because Greenland is somewhat narrower in the south, the elevations at DYE3 are 2,480 vs. 1,240 m for observed vs. modeled surface elevation, respectively. Thus, while all modeled sites are somewhat too low, the discrepancy is most significant in the south at DYE3. This may account for some, but not all, of the larger sensitivity of δ 18 O at DYE3 during DO events.
The suite of simulations shows a range of δ 18 O values and sea ice states ( Fig. 2A and SI Appendix, Fig. S1). At NGRIP, the stadial climate δ 18 O varies from −37.5 to −29‰. Of 32 initial simulations, 15 exhibit DO-type abrupt warming over Greenland with δ 18 O increases of 2.0‰ or more, where a 2.0‰ threshold could be considered representative of the minimum for a small DO-type abrupt warming (45). Each simulation with a significant δ 18 O jump starts from δ 18 O values at NGRIP of −32‰ or below (SI Appendix, Fig. S1). Most of the larger δ 18 O jumps occur under an interstadial-type salt flux (i.e., a return of salt to the North Atlantic) of size equivalent to 0.25 Sv or larger. This magnitude of salt oscillation is also in agreement with the size of salt flux required to induce significant AMOC variations in HadCM3 under present and future forcing scenarios (37). Note that there is a mean offset in Greenland between the mean model and mean data in δ 18 O of +8‰ (Fig. 1 C-E). This offset is in the same direction as in earlier isotopic simulations alongside a similar warm model bias in the temperature and a wet bias in the precipitation (21).
Temperature and Sea Ice Coefficients. There are two main approaches that are used to calculate δ 18 O-temperature coefficients and similarly, δ 18 O-sea ice coefficients. Note that the terms gradients and slopes are also used to denote these coefficient values. The first and most common approach to calculating the coefficients is to fit a second-degree polynomial to a set of stadial and interstadial δ 18 O and temperature values, minimizing the least squares term) (46,47). Coefficients (or gradients or slopes) from this method are shown in Fig. 2 C and D (bold black crosses) and SI Appendix, Figs. S2 and S3.
The second approach is to calculate the coefficients for each individual simulated DO event so that the coefficient in this case is simply ∆δ 18 O/∆ temperature for the paleothermometer. Additionally, ∆δ 18 O/∆ Arctic sea ice area indicates the sea ice coefficient. Note that Arctic sea ice area is a total Northern Hemisphere value that is calculated by summing each grid box multiplied by the sea ice concentration and grid box area. This approach yields a distribution of coefficients and is arguably a more accurate depiction of the relationship between ∆δ 18 O, ∆ temperature, and ∆ Arctic sea ice area for the suite of simulated DO events. This approach is used to characterize the distribution of paleothermometer and sea ice coefficients (Fig. 2 B-D, colored and boxplot  results).
Using the second approach and the same set of the largest DO events, we obtain a mean paleothermometer coefficient of 0.63‰ per Kelvin (16th to 84th percentile is 0.53 to 0.77) at DYE3; GISP2 yields 0.31‰ per Kelvin (0.20 to 0.46), GRIP is 0.30‰ per Kelvin (0.19 to 0.45), NGRIP is 0.29‰ per Kelvin (0.19 to 0.52), and NEEM is 0.23‰ per Kelvin (0.08 to 0.52). If the first approach is used, a similar geographical pattern and similar values emerge (Fig. 2C). This same pattern emerges independent of whether all or subsets of simulations are used in the calculations, although eliminating the smaller DO events does tend to raise the average values of coefficients. While the average simulated NGRIP coefficient is lower than the overall coefficient of 0.52‰ per Kelvin suggested for NGRIP (19), other measurements imply a somewhat lower coefficient (12), and our set of simulations indicates that the observed coefficient of 0.52 at NGRIP occurs around the 84th percentile value (Fig. 2C) (i.e., within the central range of values of simulated coefficients). Our range also encompasses the interevent set of approximate initial NGRIP coefficients from 0.28 to 0.42‰ per Kelvin used by ref. 19. Table Values. The data in Table 1 are mean values from the subset of simulations that have significant (2.0‰) DO increases in δ 18 O at NGRIP. The uncertainties are ±1 SD calculated from this sample size of 15 DO events.

Uncertainties on
Isolating the Impact of Changes in Seasonality. To qualify the relative impact of the precipitation vs. δ 18 O seasonal changes, we isolate the impact of changes in δ 18 O due to changes in the seasonal cycle of precipitation and seasonal changes in δ 18 O (34) (Fig. 3 D-F and SI Appendix, Fig. S4 where the summations are over 12 mo (index j). The superscript stadial indicates values from the cool preceding stadial-type simulated climate, and no superscript indicates values from the postwarming interstadial-type climate. This method is somewhat different from that used to decompose δ 18 O changes in ref. 48, where daily outputs were used. In particular, residuals [i.e., R esidual = (∆Pseas + ∆δ 18 Oseas) − ∆δ 18 O] could indicate changes in the higher (than monthly)-frequency covariance between precipitation and δ 18 O not captured by this monthly "seasonal"-type decomposition. However, the R esidual is smaller than ±2‰ over Greenland (less than ±0.5‰ over most of central Greenland). The stadial to interstadial values are calculated using 50 y of data from each climate, with the 20 y on either side of the abrupt warming excluded from the analysis.