Decadal-scale progression of Dansgaard-Oeschger warming events

During the last glacial period, proxy records throughout the Northern Hemisphere document a succession of rapid millennial-scale warming events, called Dansgaard-Oeschger (DO) events. A range of different mechanisms have been proposed that can produce similar warming in model experiments, however the progression and ultimate trigger of the events is still unknown. Because of their fast nature, the progression is challenging to reconstruct from paleoclimate data due to the limited temporal resolution achievable in many archives and cross-dating uncertainties between records. Here we use new 5 high-resolution multi-proxy records of sea-salt (derived from sea spray and sea ice over the North Atlantic) and terrestrial (derived from the Central Asian deserts) aerosol concentrations over the period 10-60 ka from the Greenland NGRIP and NEEM ice cores in conjunction with local precipitation and temperature proxies from one of the cores to investigate the progression of environmental changes at the onset of the warming events at annual to multi-annual resolution. Our results show on average a small lead of the changes in both local precipitation and terrestrial dust aerosol concentrations over the change in sea-salt 10 aerosol concentrations and local temperature of approximately one decade. This suggests that, connected to the reinvigoration of the Atlantic Meridional Overturning Circulation and the warming in the North Atlantic, both synoptic and hemispheric atmospheric circulation change at the onset of the DO warming, affecting both the moisture transport to Greenland and the Asian monsoon systems. Taken at face value, this suggests that a collapse of the sea-ice cover was not the initial trigger for the DO events. 15

These records are of great value to the scientific community, and the analysis is meaningful and appears to be rigorously done (request for minor clarification below).This paper is clearly a great contribution to the literature, and I have only some minor suggestions that may improve the clarity and interpretation.

Thank you!
(1) The current manuscript only describes the relative phasing of the onset, midpoint and endpoint of each transition.What is missing is an analysis of the duration of each of the transitions, to put the lead/lag values into perspective.For example, if the transition were to take 100 years, then a 10-yr lead indicates that the climatic components reflected by these records co-evolved; if the transition were to take 5 years only then the same 10-yr lead suggests a decoupling, with the shift in one component (say the jet) completed before the others respond.I request the authors add one panel to Figures 3 and 4 each that gives the transition duration.
We agree that this is an important point.Because both Figure 3 and 4 show difference relative to Na, we prefer not to add absolute information to these figures.However, to answer the reviewer's comment, we have added the relevant information to the text (P 7, L 11 ff): "For all of the transitions, the inferred timing differences relative to the onset of the transition in Na are smaller than the duration of the transition itself in each of the proxy records.That means that none of the proxies exhibit a complete stadial-interstadial transition before the onset of the transition in the sea-salt aerosol concentration."And later in the discussion (P 12, L 25 ff): "The fact that for all transitions the inferred timing difference relative to the transition in Na is smaller than the duration of the transition in that parameter indicates that the respective parts of the climate system co-evolved over the transition.That means that the changes in atmospheric circulation at the DO-onset where not completely decoupled from the change in sea-ice cover and Greenland temperature."Many of us are visually oriented.Would it be possible to show the D-O average transition in Ca, Na, lambda and d18O together in one plot like is done in Fig. 2 for Ca? Either the data or just the fits (if the data are too messy).This would really give a nice visual representation of how the "average" transition occurs.Because all of the transitions have different lengths it is not possible to meaningfully stack them.This is the reason why we chose to combine the estimates of the relative onset timings and not the complete transitions.We would like to point out that Figure 2 only shows the fit of a single transition, not a stack.
(2) The paper only analyzes the glacial-to-interglacial (or D-O warming) transitions, and not the interglacial-to-glacial (or D-O cooling) transitions.Have you tried analyzing the latter?I would be very interested to see the phasing for these transitions also.I imagine it is more challenging given the smaller and less abrupt nature of these transitions but I think it would be valuable nevertheless.
We agree with the reviewer that analyzing the cooling transitions would be very valuable.However, in the presented manuscript, we deliberately chose to focus on the DO onsets as they are much more clearly distinguishable in the records and are possible to describe using the relatively simple model employed here.Furthermore, the fundamental difference between the warming and cooling transitions in their manifestation in the records suggests a different mechanism and thus merits its own separate investigation which we plan to do in the future.
(3) How precise is the depth registration of the various CFA components relative to one another?And relative to the layer thickness and d18O records?I can imagine there may be cm-scale offsets, which could become important given the extremely small time phasings that the authors interpret.Please address this briefly.
Between the different CFA records, the depth assignment is accurate to within a few mm as they are measured on the same piece of ice and the relative component alignment is measured regularly as described in the cited references.Because the NGRIP CFA records form the data basis for the annual layer-counting, the CFA depth scale is the fundamental depth scale to which the layer thicknesses and dating refers.The relative uncertainty of the depth assignment of the CFA/layer thickness data to the d18O record is fundamentally limited by the accuracy of the subsampling within each 55 cm piece of ice (called a bag).Typically, sample depths are measured with a meter stick and are accurate to within half a centimeter or better.We added the relevant information to the manuscript.(P 3 L 29 ff): "Between the individual data sets the co-registration uncertainty is limited by the absolute depth assignment of the data sets.This uncertainty is typically on the order of a few millimeter to at most single centimeters and thus on the which translates to a co-registration uncertainty in the sub-annual range."(4) Appendix A was not completely clear to me, and I think it should be elaborated on in some more detail for the lay reader.Do I understand correctly that the function that is being fit to the data is the ramp function plus an AR(1) noise function, and that the parameters of both are varied in the Monte Carlo sampler?
The language of prior and posterior distributions suggests this is a Bayesian approachplease confirm.
How is the goodness of fit evaluated, and what criterion is used in the MC sampler to accept or reject individual solutions?It would help greatly if you could add a figure showing how several individual iterations of the fitting process look like.
We believe that a thorough introduction to Bayesian statistics and Markov Chain Monte Carlo methods is beyond the scope of this paper.Nevertheless, we want to give an intuition on the sampling process here: In the sampling process, the sampler randomly walks the whole parameter space, in this case represented by the four parameters describing each transition and the two noise parameters.At each step, the sampler evaluates the likelihood of observing the data given the current set of parameters and accepts or rejects these parameters with a probability proportional to this likelihood.If performed over a large number of samples, this process will yield samples from the posterior distribution which is what we use for our analysis (see for example Gelman et al., 2013) Gelman, A. et al. (2013), Bayesian data analysis, Chapman and Hall/CRC (5) I am confused by the statements below Fig 4 (Starting on P9 L18).If both Ca and lambda lead Na by ∼10 years, how come these two are not necessarily synchronous?This is very counter-intuitive; all records are evaluated on the same depth scale, so why wouldn't they be?I think the relative phasing of all these records is the central result of this paper, so it would be important to establish a robust sequence of events.What would be needed to establish synchroneity of Ca and lambda?Would you need to run the analysis again relative to the Ca transition instead of relative to the Na transition?If not too much work, that may be worth doing, given the importance for the interpretation.
I could imagine a 4x4 matrix for NGRIP with the lead/lag of each of the records evaluated relative to the others, and a 2x2 matrix of the same for NEEM.
We agree that this is indeed counter-intuitive.We slightly rephrased the sentence to make this point clearer (see below) The reason for this is that the uncertainties both for the lags in the individual transitions as well as for the combined evidence is very strongly dependent on the estimates from the layer thickness timeseries.Depending on how these highly correlated uncertainties are projected onto the different axis (lt-Na, Ca-Na, Ca-lt) slightly different distributions arise.For Ca and lt, the timing difference at the onset of the transition shows a lead of Lambda over Ca with 4 (-5/+4) yr that is not incompatible with a zero lag at the 95% probability level.A different way of looking at the order of events is to calculate the average probability of the order of the transition onsets.This shows that the transition Lambda and Ca are about equally likely to occur first, whereas Na and d18O are about equally likely to occur last.From this type of analysis, one can also establish the most likely order of onsets during the warming transitions which is lt, Ca, d18O, Na for NGRIP and Ca, Na for NEEM.We have added these results in a new Figure focusing on the NGRIP ice core into the revised manuscript (P9, L16 ff and Figure 5)."Note that the density functions shown in Figure 4 cannot be used to infer timing differences between the other parameters.This is a direct result of the estimates being conditional on the timing of the transition in sodium, leading to large correlations between the lag estimates for the other parameters.That means that even though e.g two probability density functions of the differences relative to the transition in sodium largely overlap, that does not necessarily mean that the timing difference is equal to zero.In the case of the timing difference between the transition onsets of the increase in annual layer thickness and the decrease in Ca 2+ concentrations the combined lead of the change in annual layer thickness is not larger than zero at the 0.95 probability level with 4(-5, +4) years.To establish the most probable sequence of events at the transition offset we calculate the average order of the onset times, shown for NGRIP in Figure 5.The average positions show that the change in accumulation and Ca 2+ concentrations about equally likely occur first whereas the transitions in Na + and δ 18 O about equally likely occur last.The same analysis for the NEEM results confirms this sequence."(6) All age axes have a "BP" label.Do you use BP 1950 or b2k?This is a contentious point in the ice core community (Wolff, 2007), but BP 1950 is the best choice in my view based on precedent in the literature and convention of the radiometric dating communities.At least specify which is used.
We had only specified this in the caption of Figure 1, but we have now also added this to the text (P4, L1) "All ages are given relative to 1950."We have also changed all axis-labels to "before 1950" instead of "BP" to avoid confusion with radiocarbon ages.

Line-by-line comments:
P1 L17-19: The phrase "DO event" is unfortunately ambiguous, with some people equating them to the abrupt warming phases, and others to the interstadials.To avoid this, consider ". . .millennial-scale abrupt climate change, called the DansgaardOeschger cycle (REFS).During abrupt DO warming, . .." Thank you for the suggestion.To make the text more consistent we changed the wording throughout the text so to specifically state that we are looking at the onset of the D-O events.However, we deliberately do not use the term "DO cycles" to avoid any notion of cyclicity in these events.L22: I assume your are still talking about DO warming here?Please specify that the changes described are for the warming phase.
Clarified in the text.P2 L6: Consider replacing Henry et al. with the recent review by (Lynch-Stieglitz, 2017), to avoid arbitrarily picking one study out of dozens that demonstrate the link to AMOC.

Done.
P2 L23: "Some of . ..".I think you can safely say "All of" (or "most of " to be conservative).I am not aware of any model simulation or theory that does not involve sea ice as either the trigger or amplifier.You simply cannot get that much warming that quickly without sea ice.
We have now changed the sentence and use "Many of…".Most of the experiments are run without an interactive sea-ice model and thus cannot realistically inform on sea-ice feedbacks.We do however agree that most of the recent model studies with coupled models involve the sea ice as triggers or amplifiers.P3 L20: "exact co-registration".How exact?Please specify relative and absolute depth registration of various CFA components.
We added the relevant information to both the CFA data sets as well as to the end of the paragraph: "…exact co-registration of the aerosol concentration records at the millimeter scale…" (P3, L21) "Between the individual data sets the co-registration uncertainty is limited by the absolute depth assignment of the data sets.This uncertainty is typically on the order of a few millimeter for CFA data up to around a centimeter between d18O and CFA data, which translates to a co-registration uncertainty of sub-annual to annual range."(P3, L29 ff) P3 L23: Do you use the actual single layer annual counts, or the 20yr averages that are publicly available?
The data presented here is based on the single annual layer counts and new 10yr averages are published alongside this paper for the whole 60 ka as well as higher resolution datasets for the individual transitions.To make all axis directions consistent, we changed the axis in Figure 2 to match those in the other figures.P5 L30-31: So are your interpreting the changes in terms of the source strength only?Or does transport to the site dominate?Some would argue for the latter.
The exact interpretation in terms of source and transport changes is described in detail in the discussion section of the paper.(P10 ff) P5 L33: Do you actually fit an exponential, or do you fit a linear to the log(Ca) time series?
We fit a linear function to the log(Ca) and log(Na) time series, as described in the appendix.We added this also to the main text for clarification: (P6, L2 ff) "…or exponential transition (i.e. a linear transition fitted to the log-transformed data for all other records) between two constant levels…" P6L16: "decreases" should be "increases" here, right?(it goes from 2 to 7yr, so increase?) We prefer decreases in this case as "high resolution" typically referrers to lower sampling intervals and vice versa.
Figure 2: please add units to the y-axis labels.Also, how come the fit is so smooth / rounded?Isn't the fit a linear ramp?Is this because you average over many solutions?
Yes, the smoothness is a feature of the marginal posterior median ramp which is shown here.The individual realizations are linear ramps with sharp kinks.We have now added units to the y-axis labels P9 L5: Buizert et al. 2015 should be cited as WAIS Divide Project Members 2015.
We have now updated the reference accordingly.P10 L23-24: Add a reference for this claim.
We have now rephrased to "…lack of synchronous changes…" P10 L32: "events" is confusing here.Are you talking about individual synoptic / precip events?or DO events, better specify more clearly.
We have now specified as "precipitation events".P11 L15: This idea was suggested by (Seager and Battisti, 2007) and (Wunsch, 2006) We have now added Wunsch, 2006 as original reference.P11 L16: I had expected a larger discussion about wet vs. dry deposition.Could the coincidence of lambda and Ca changes be explained that way to some degree?Ca is, as sea salt, very efficiently scavenged by snowfall events and its deposition is hence governed by the frequency of the precipitation events.P11 L18: effect should be affect Corrected, thanks for pointing this out.
P11 L34: This further supportS. . .Corrected, thanks for pointing this out.P12 L13: "reduction of the sea ice cover that ultimately coincided with the Greenland warming AND WAS PRESUMABLE A MAJOR DRIVER THEREOF" Again, I think it's hard (impossible?) to get such a large Greenland temp response without a change in sea ice cover.
We completely agree and we do not argue against the fact that the sea-ice change is likely the major driver of the warming in Greenland.P18 L25: What is the rationale for taking the log of lambda instead of just lambda itself?
We have now added the missing symbols to the text.

Reply to reviewer comment #2
In this study the authors use high-resolution ice-core records of aerosols, water isotopes, and layer thickness from Greenland to examine phasing of different aspects of the climate system during Dansgaard-Oeschger events.They use objective techniques to tease out small leads and lags between the noisy records, showing that changes in calcium aerosols and layer thickness lead changes in sodium aerosols and water isotope ratios.They conclude that these lags suggest that changes in sea ice extent did not occur before other changes in the climate system, namely atmospheric circulation.The manuscript is very well written, the analysis is careful, and the discussion is well-argued.The results are quite impressive and should be of wide interest to the community.I have a couple questions and concerns that I hope the authors will address and then several minor questions.
Thank you! (somewhat) Major questions: The authors have made a compelling observation through their analysis of leads and lags between the proxy records.They attribute these differences in timing to aspects of the climate system that have different influence on the proxies.My first two questions have to do with whether these leads and lags could arise from other factors.I suspect that these concerns are not too important to the conclusions of this study.
1. How does the signal-to-noise ratio of the record influence fitting the transition model and the identification of starting, mid, and end points?Here I'm thinking of the SNR as quantified by the size of the transitions compared to the variance within the stadials and interstadials.I realize the fitting procedure takes into account the interannual noise and its autocorrelation.But if you have an idealized known ramp function with different levels of background noise, will the model find the same starting, mid, and endpoints?One could imagine that an increased background noise could lead to the identification of time-shifted transition points depending on the fitting technique.One then risks conflating difference in the timing of signals between records with difference in one's ability to detect signals between records.I cannot tell from the description of the transition model alone how much of an issue this is to this analysis.Doing some simple tests with a couple different ramp-fitting and significant change detection techniques (though ones less sophisticated than the technique used by the authors), I find that different levels of interannual noise can influence the timing of a fitted transition, though not in all circumstances.
My particular worry here is that the substantially higher noise (interannual variability) in the Na records could lead to the identification of a delayed onset or shorter transitions compared to the Ca and other records.My worry is heightened slightly in that the relative timing seems to correspond to the level of background noise (at least visually): the d18O and Na timing are most similar among the proxies and also both appear to have much lower SNR (more noise) compared to the Ca and layer thickness records.Do the mean lags depend on the amplitude of the DO event or the length of the transition between stadial and interstadial?This could be the case if the ramp-fitting depends on the amplitude of back ground noise.I realize this could be hard to determine since one needs to look at many events at once to see the mean lags.But a scatter plot of lags vs. event magnitude or fitted ramp duration could be informative.I suspect that this concern is entirely accounted for by the very careful analysis of the marginal posterior densities of the onsets, midpoints, and endpoints for each proxy and the comparison between proxies, that the authors have already performed.It would however be helpful to see the influence, or the demonstration of the absence of influence, of the SNR on the fitting procedure given that the conclusions rest on the difference in timing with respect to Na in particular.I'd find it very useful to see this demonstrated on artificial ramp signals (of varying duration) where the true timings are known explicitly, with varying SNR, and especially with the SNRs relevant to d18O, lambda, Ca, and Na.It seems crucial to know that different SNR alone can not account for the difference in timing identified by the fitting procedure.
We do completely agree with this concern and thank you for bringing this up.You are correct in that the signal to noise ratio as defined by the amplitude of the ramp divided by residual standard deviation varies widely between the different parameters and that Ca has by far the highest SNR of the records.In our analysis this is reflected in the posterior standard deviations of the timing estimates for the individual parameters: The higher the SNR, i.e. the clearer the transition, the lower the posterior standard deviation.To illustrate this, the Figure below shows the marginal posterior standard deviations of the start point of the transitions as functions of the SNR for the NGRIP record.
As we propagate these uncertainties to the lead/lag calculations and the combined estimates we take this into account.This is also clear when plotting the leads/lags relative to Na as a function of the SNR in the respective parameter which does not show any relationship with the SNR: We have added this discussion and the two Figures above to the Appendix of the manuscript.
2. Can the authors rule out the influence of water isotope diffusion on the difference in timing of the d18O signals and those of Ca and layer thickness?Based on analysis of NGRIP (Gkinis et al 2014) and a similar site in Antarctica (Jones et al 2017(Jones et al , 2018)), I'd guess diffusion lengths are on the order of 5-10cm through this interval and so are not insignificant compared to the annual layer thickness.Such diffusion lengths can have meaningful influence on the inter-annual and even decadal variability in the water isotope record (Jones et al 2018).I imagine that correcting the records used here for the potential influence of diffusion, if that would even be sensible, is far beyond the scope of this study.However, it seems entirely reasonable to estimate the influence (if any) of the smoothing implied by diffusion on the timing identified by the transition model fitting procedure.If you take identical idealized ramps, and smooth one with a time-scale reflective of water isotope diffusion lengths, will the fitting procedure identify the same start, mid, and end points?I suspect these effects, if any, are small, though we are only talking about lags of a few years.
Yes, the smoothing by the isotope diffusion would influence the shape of the transition.It is however difficult to say what the exact effect on the estimated onsets would be as it also increases the SNR by reducing the high-frequency variability that is captured by the noise term.That being said, in the case of DO onsets, where the d18O signal rapidly increases, any diffusion/smoothing of that signal would yield to a slight shift towards earlier onsets, i.e. would reduce the lag between Ca and d18O.Nevertheless, this is an interesting and important thing to check.To do so, we examined the auto-correlation time that is estimated by the fitting algorithm.Any diffusion of the signal would lead higher auto-correlation times estimated for d18O than for all the other parameters, even though the AR(1) process is not an ideal description for the autocorrelation introduced through gaussian smoothing.
This Figure shows these estimates for all parameters in NGRIP.In comparison to Ca, d18O exhibits similar auto-correlation lengths throughout the record.Indeed all the records show an increase in auto-correlation length with increasing depth with is a direct consequence of the decrease in the resolution in terms of years due to the layer thinning.
We have added the following to the discussion (P12, L29 ff) "Due to the small timing differences between the records it is worth noting that water isotope records from polar ice cores are subject to smoothing by diffusion.For the NGRIP isotope record the diffusion length at the end of the last glacial is on the order of five to ten centimeters (Gkinis et al., 2014), influencing the high frequency variability.For the analysis here the diffusion means that the rapid increase of the δ 18 O signal at the D-O onsets would be slightly shifted towards earlier times, leading to lower apparent lags between the aerosol records and the water isotope record.Thus, the inferred lead of Ca over d18O can be regarded as a conservative estimate."Gkinis, V. et al. (2014), Water isotope diffusion rates from the NorthGRIP ice core for the last 16,000 years -Glaciological and paleoclimatic implications, Earth andPlanetary Science Letters, 405132-141, doi:10.1016/j.epsl.2014.08.022 3. The most interesting conclusion of this study to my mind is the authors statement that "at face value, this sequence of events suggests that the collapse of North Atlantic sea-ice cover is not the initial trigger for the DO events..." Because of its potential wide interest, this statement deserves some scrutiny.It rests on the authors use of Na as a "qualitative indicator of sea ice cover" in the North Atlantic.
The discussion on Page 5, lines 3-28, highlights the debate over in the interpretation of Na very well.The authors interpretation is laid out on Page 5 lines 20-25.Markle et al 2018 find that most of the millennial variability in Antarctic sea salts can be explained simply by the changes in moisture rainout that are required to explain the water isotope record (these changes also explain most of the changes in Antarctic Ca variability and its relationship to water isotopes in both Antarctica and Greenland (c.f.their Figure 4)).This suggests comparatively small if any changes in sea salt source latitude or strength are needed to explain those Na records (though changes in those things are still possible of course).Is there evidence that this explanation for the sea salts is insufficient in Greenland?Are the observed changes in sea salt for example much larger than what one would expect from temperature dependent rainout alone?It was unclear to me from this discussion that the sea salt source strength (or even mean source latitude) should have a clear relationship to the sea ice edge.
Generally, for Greenland the same mechanism as described by Markle et al. (2018) explains most of the millennial-scale variability as shown in Schüpbach et al. (2018).However, applying the model of Markle et al. (2018) to the DO events would completely preclude any leads or lags between d18O and the aerosols as well between the different aerosols.On shorter timescales dynamical changes in the (co-) transport of aerosols and moisture seem to play a relatively larger role in determining the concentrations observed in the ice cores.Here, this results in the lack of co-evolution between the investigated proxies and can explain the lack of coherence between the records in Markle et al. 2018 at the sub-millennial timescales.Schüpbach, S. et al. (2018), Greenland records of aerosol source and atmospheric lifetime changes from the Eemian to the Holocene, Nature Communications, 9, 1-10, doi:10.1038/s41467-018-03924-3.
If the main way sea ice influences Greenland Na is through its relationship to the variables driving the rainout effect, then a change in sea ice extent doesn't necessarily mean one should have a coincident change in Greenland Na, particularly at interannual timescales.
For example, one can imagine a scenario in which sea ice extent begins to retreat at exactly the same time as the changes observed in Ca initiate.Coincident increases in temperature and moisture removal would cause a decrease in the amount of Na (and Ca) reaching Greenland (as described by the authors).However, if that change in sea ice extent caused an increase in sea salt source production or a northward migration of mean source latitude (both debatable but reasonable, particularly the latter) this could temporally compensate for the increased removal.This combination of influences could lead to an apparent timing difference in the final Na signal observed in Greenland compared to the actual timing of sea ice changes (an example of the superposition of competing source and rainout factors on polar aerosols is given in the Supplement of Markle et al 2018 c.f. Figures S9).A somewhat similar scenario may be likely for the water isotopes, as changing sea ice extent could drive moisture source effects that could temporarily compensate for the decreased depletion driven by simultaneous changes over the ice sheet (these would be of the correct sign to compensate, though it would be hard to assess the potential size of this effect on the isotopes without analyzing deuterium excess records from the same cores).Quantitatively disentangling these influences may be well outside the scope of this study.But this does at least suggest a limitation to using Na as a qualitative indicator of sea ice, and should suggest some commensurate caution in the conclusions drawn based on that interpretation.
Even in the absence of competing influences, uncertainty in the linearity of the sodiumas-seaice-extent interpretation poses challenges.Dose a change in the sea ice edge (or extent) of a given size have the same impact on Na in Greenland if the sea ice edge is at 55 degrees North (just for example) versus if the edge is at 65 degrees North?If relationship between changes Greenland Na and the sea ice absolute position is (sufficiently) nonlinear then the changes in sea ice at the start of DO events may not be as detectable in Greenland Na as changes later in the event.This could lead to apparent lags of the Na signal to Ca, even if the change in sea ice initiates at the same time as the change in whatever drives Ca.I'm not sure I think this is particularly likely, but it is certainly plausible.Further, the impact of sea ice on climate isn't linear.Again a given size change in the sea ice edge when the edge is at 55 degrees North (for example) doesn't have the same impact on the surface radiation budget nor the buoyancy forcing of the overturning circulation as when the edge is at 65 degrees North.The relative timing of sodium changes in Greenland don't necessarily rule out sea ice changes as a potential "trigger" of DO events, if climatically meaningful sea ice changes can happen without influencing Greenland Na.
To be sure, there are innumerable, potentially ad hoc, explanations for the data, and it is not the authors responsibility to come up with and then weigh the merits of all such explanations.However the assertion that the lag in sodium seen in Greenland necessitates a lag in changes in sea ice extent with respect to shifts in atmospheric circulation, and that this in turn rules out sea ice changes as the trigger for DO events, is somewhat bold.That this "provides an essential benchmark for climate models" is bolder still.Both of these statements to my mind somewhat outpace the robustness in the interpretation of Na as a qualitative indicator of sea ice extent.I think some slight tempering of the conclusions is merited here and perhaps some discussion of the limitations to the interpretation, or much more compelling evidence is needed that a change in sea ice extent must be seen in the Greenland Na records regardless of other processes.To be clear, I do think the difference in timing of the aerosol species identified here is a compelling target for modeling.And I do think the climatic interpretation offered by the authors is plausible, and one to be taken seriously.It is just not clear to me that these data alone place that strong of constraints on the timing of sea ice changes.
We agree with the possible complicating factors you pointed out.However, it is difficult to imagine that a "climatically meaningful" sea-ice change can happen without influencing neither Na nor d18O in Greenland ice through the changing source signals, en-route rainout or changes in seasonality or a combination of these factors: Taking the transient DO simulation performed by Vettoretti andPeltier (2015, 2018) as an example, the change in the sea-ice coverage in the North Atlantic results in an immediate increase in moisture cycling (evaporation and precipitation) over the open ocean as well as an increase in temperature in Greenland.Nevertheless, we agree that barring any direct simulations using both isotope and aerosol models over DO-events, the direct test of our hypothesis will be very difficult.To account for this limitation, we have rephrased the relevant statements in the conclusions to be more careful (P13, L11 ff) "The progression of environmental changes revealed in the Greenland aerosol records provides a good target for climate models that aim to transiently simulate DO events, preferentially explicitly modeling both water isotope and aerosol transport."Vettoretti, G., and Peltier, W. R. (2015), Interhemispheric air temperature phase relationships in the nonlinear Dansgaard-Oeschger oscillation, Geophysical Research Letters, 42,1180-1189, doi:10.1002/2014GL06289 Vettoretti, G., and Peltier, W. R. (2018), Fast Physics and Slow Physics in the Nonlinear Dansgaard--Oeschger Relaxation Oscillation, Journal of Climate, 31, 3423-3449, doi:10.1175/JCLI-D-17-0559.1 Minor points: Page 1, line 8: The authors say "from one of the cores".... Which one?
We have now clarified the sentence in the abstract.
Page 1 line 17: I think that the clause "In the course of the last glacial period" should be moved after the 2nd comma in this line (after "warming episodes,").This is a possibly silly language thing but: the ice core records reveal that the events were in the last glacial period.It wasn't during the last glacial period that the ice core records revealed these events (that was in the 1980s!).

Done.
Page 1 line 20: This is a language choice the authors should feel free to ignore but "..going along with an almost doubling..." sounds strange to my ear, though I can't tell why."...coincident with a near-doubling..." sounds better to me.

Adjusted accordingly.
Page 2 line 3: I'd move the "Also" at the start of the sentence to after "Northern Hemisphere".

Done.
Page 2 Lines 20-35: This is a very nice description of several related though distinct ideas.Thanks!Thank you!Page 3 lines 15-20: Are the Ca and Na records corrected for the sea-salt vs non sea-salt components of each (e.g. using average Na/Ca mass ratios in average crust and sea salt, as is common in Antarctic records)?I ask because having looked into this once briefly, it seemed like the corrections often used in Antarctica lead to nonsensical results in Greenland, which was disturbing, and I wasn't sure why.
No, the records we use here are not "corrected" for the sea-salt and dust contributions.We chose not to do this for two reasons: First and foremost, the Ca/Na ratio in the dust in Greenland is quite different from the crustal average that is typically employed for Antarctic records.The reason for this is the large Ca content of Central Asian Dust.This is also true for Antarctica and has previously been extensively investigated by Bigler et al (2006).Page 6, Line 9: I think there is a missing "of" between "interpretation" and "phase".
Page 6, lines 29: If you don't make the assumption that the DO-events have the same imprint in both cores, how would that effect any of your conclusions?
Making this assumption would allow to combine the estimates from the two cores into one.We chose not to do this to be able to gauge the robustness of our estimates, i.e. whether they agree between the cores.
Page 10, lines 25-30: Might be worth citing Fudge et al 2016 here, they show a strong divergence of accumulation rate and d18O in ice cores at millennial time scales.
Thank you for pointing out this reference.We added it to the manuscript with the qualification that it deals with Antarctica and not Greenland.
Page 11, lines 3-5: This is nice, convincing analysis.To explain the complete amplitude of the stadial/interstadial changes in the Ca concentration in the ice, a source strength change of a factor around 4 is needed.We have added this to the manuscript (P12, L9) and refer the Reviewer to Schüpbach et al. (2018) for further details:
Also other ::::: Other proxy records throughout the Northern Hemisphere ::: also : document the widespread environmental imprint of these rapid warming events.Marine sediment cores (Dokken et al., 2013) and aerosol records from Greenland (Spolaor et al., 2016) show a reduction in perennial sea-ice cover in the Nordic Sea and Arctic Basin and ocean circulation proxies indicate an increase in the Atlantic meridional overturning circulation (AMOC) (Lynch-Stieglitz, 2017).Records from North America indicate a change in the moisture advection from the Pacific with dryer conditions during the warm interstadial periods likely related to changes in atmospheric circulation (Wagner et al., 2010;Asmerom et al., 2010)resulting in : .::::: These ::::::::: circulation ::::::: changes ::::::: coincide :::: with increased wild-fire activity in North America as clearly imprinted in the Greenland ice-core record (Fischer et al., 2015).Furthermore, records from Eurasia indicate rapid changes in the local ecosystems (Rousseau et al., 2017).In the lower latitudes, : speleothem and sediment records from both South America (Wang et al., 2004;Deplazes et al., 2013) and Eastern Asia (Wang et al., 2008) indicate a northward displacement of the Inter Tropical Convergence Zone (ITCZ) at the time of DO-warming ::: DO :::::::: warming : (Cheng et al., 2012) resulting in rapid changes in tropical hydro-climate and methane emissions from the tropical wetlands (Baumgartner et al., 2014).A recent synchronization of cosmogenic radionuclide records from ice cores and low-latidude :::::::::: low-latitude speleothems back to 45 ka ago shows that atmospheric circulation changes in the tropics occurred synchronously with the Greenland warming within the cross-dating uncertainties of around 180 yr (Adolphi et al., 2018).The atmospheric circulation changes associated with the Asian monsoon systems documented in the Asian speleothems also reduced the mobilization and export of mineral dust aerosol from the Central Asian deserts during the warm stadial periods as documented by downstream sediment records (Porter and Zhisheng, 1995;Jacobel et al., 2017).
A range of different, non-exclusive mechanisms that produce interstadial-like warming events during preindustrial and glacial climate conditions have been proposed and tested in model experiments.These include the direct modulation of the AMOC by freshwater addition to the North Atlantic (e.g.Knutti et al., 2004), where ceasing artificial freshwater forcing in the North Atlantic results in an increasing strength of the AMOC and subsequently increasing Greenland temperatures.Some ::::: Many of the proposed mechanisms involve the reduction of North Atlantic sea-ice cover either as driving or as amplifying process for the warming.In both coupled and uncoupled model experiments, the removal of winter sea-ice cover in the North Atlantic and Nordic Seas alone generates an increase in Greenland temperatures and snow accumulation rates similar to observations by exposing the relatively warm underlying ocean (Li et al., 2010(Li et al., , 2005)).Experiments with coupled atmosphere and ocean show that this reduction in sea-ice cover can be induced by changes in the wind stress over the sea ice either arising spontaneously (Kleppin et al., 2015) or controlled by elevation changes of the Laurentide ice sheet (Zhang et al., 2014).Furthermore, heat accumulation under the sea ice during stadials can destabilize the otherwise strongly stratified water column, eventually leading to breakdown of the stratification and melting of the sea ice from below (Dokken et al., 2013;Jensen et al., 2016).Alternatively, Peltier and Vettoretti (2014) and Vettoretti and Peltier (2015) describe the occurrence of spontaneous DOlike oscillations in their fully coupled model caused by the buildup and break down of a meridional salinity gradient between the open and the sea-ice-covered North Atlantic.The collapse of the salinity gradient then leads to a rapid disintegration of the sea-ice cover and an invigoration of the AMOC.In summary, it is clear that sea ice is a crucial factor in generating the full climatic change seen at the DO-warming ::: DO :::::::: warming (Li and Born, 2019).However, the sequence of events, and whether the sea ice ::::: sea-ice : loss is triggered by atmospheric or oceanic changes, differs between the proposed mechanisms.
In any case, the progression and possible interaction between the distinct processes that can partake in a DO event :: the :::: DO ::::::: warming : are difficult to constrain from paleo-observations and their validity therefore difficult to test, especially as some of the proposed mechanisms lead to virtually indistinguishable model results (Brown and Galbraith, 2016).The relative timing between the start of DO events in proxy records of different parts of the Earth system can yield critical insight into their spatiotemporal progression and possible causal relations.However, due to the fast onset of the DO events, relative time differences are expected to be in the order of years to decades, requiring very high temporal resolution of paleo-records and often unattainable small relative dating uncertainties between records from different archives.
For both ice cores, the current versions of the GICC05 age scale were used, which in the case of NEEM has been transferred from NGRIP using volcanic match points (Rasmussen et al., 2006;Andersen et al., 2006;Svensson et al., 2006;Rasmussen et al., 2013).:: All :::: ages ::: are ::::: given ::::::: relative :: to ::::: 1950.Even though at the volcanic match points the dating uncertainty is relatively small, the uncertainty introduced by the interpolation between the match points generally precludes direct comparison of absolute timing between the two cores with the precision required for this study.O records (f).All records are shown in :: as decadal averages on their respective version of the GICC05 age scale relative to 1950 (Svensson et al., 2008;Rasmussen et al., 2013).The vertical lines mark the investigated interstadial onsets as given by Rasmussen et al. (2014).Note that the vertical scales for the aerosols and annual layer thickness are logarithmic.
Mineralogy and isotopic composition of mineral dust aerosol from central Greenland indicate that its main :::::::: dominant :::::: glacial sources are the central Asian Taklamakan and Gobi deserts (Biscaye et al., 1997;Svensson et al., 2000).Large dust storms occurring during spring can lift dust up into the westerly jet stream where it is transported over long ranges and at high altitudes (Sun et al., 2001;Roe, 2009).The dust emission and subsequent entrainment into the jet is strongly dependent on the specific synoptic circulation in the area which itself is governed by the latitudinal position of the westerly jet (Nagashima et al., 2011;Roe, 2009).Under current conditions, the position of the westerly jet over central Asia varies seasonally, changing from south to north of the Tibetan plateau during late spring, in unison with the northward movement of the East Asian Summer Monsoon rain belt (Schiemann et al., 2009;Yihui and Chan, 2005).During the last glacial, the jet was located further south, and, especially during cold periods, enabling more frequent or even permanent conditions for the deflation and entrainment of the central Asian dust (Chiang et al., 2015;Nagashima et al., 2011).Once in the jet stream, the mineral dust aerosol is transported above the cloud level and is largely protected from scavenging by precipitation and transported efficiently to Greenland (Schüpbach et al., 2018).This allows us to interpret changes in calcium concentration records in terms of changes in the conditions needed for dust entrainment, i.e. the latitudinal position of the westerly jet and the Asian hydro-climate.
Sodium-containing sea-salt aerosols are produced either by bubble bursting at the surface of the open ocean or blowing saline snow on the surface of sea-ice (Wagenbach et al., 1998) ::::::::::::::::::::::::::::::::::: (Wagenbach et al., 1998;Yang et al., 2008).Sea-salt aerosol is transported together with moist air masses from the North Atlantic onto the Greenland ice sheet (Hutterli et al., 2006).Under current conditions, aerosol transport models point at the emission from blowing snow from winter sea ice as the dominant source of sea-salt aerosol for Greenland in winter and suggest that this source is responsible for a large fraction of the seasonal variability seen in ice-core records (Huang and Jaeglé, 2017;Rhodes et al., 2017).However, on inter-annual time scales the influence of the atmospheric transport and deposition dominates the variability at the central Greenland sites for recent times, due to the overall low contribution of sea-ice-derived aerosol to the total sea-salt aerosol budget (Rhodes et al., 2018).
Under glacial conditions, the extended multi-year sea-ice cover moves both sources of open-ocean and sea-ice-derived seasalt aerosols further away from the ice-core sites.Intuitively, this would lead to a reduction in sodium deposition on the ice sheet due to the longer transport.However, the opposite is observed in the Greenland ice-core records of the last glacial with much higher concentrations during cold climate periods than during warmer periods (Schüpbach et al., 2018;Fischer et al., 2007).Isolating the effect of the more distal sources : , Levine et al. (2014) have shown , that the more extensive sea ice around Antarctica under glacial conditions would lead to reduced sea-salt aerosol transport to the ice sheet.This effecthowever : , ::::::: however, : is overcompensated by changes in the atmospheric circulation that enhance production and transport during the glacial in comparison to present day conditions, leading to observed increase in sea-salt aerosol concentrations in Antarctic ice cores (Wolff et al., 2010).This suggests that changes in the transport and deposition regimes must be responsible for a large fraction of the glacial/interglacial and stadial/interstadial variability observed in the Greenland ice-core records (Schüpbach et al., 2018;Fischer et al., 2007).One plausible mechanism to explain ::::::::: explanation ::: for : the apparently more efficient transport of sea-salt aerosol to Greenland during the cold periods are the drier conditions over the cold, sea-ice-covered North Atlantic compared to open-ocean conditions.This implies that changes in the sea-ice cover do affect the sea-salt aerosol concentrations in the Greenland ice cores not only because of an influence on the source of the sea-salt aerosol but also on the efficiency of its transport, i.e. its deposition en-route.Reduced sea-ice cover increases evaporation from the open ocean, resulting in increased scavenging en-route and subsequently reduced transport efficiency of sea-salt aerosol to the Greenland ice sheet, especially because of its co-transport with moist air masses from the North Atlantic (Hutterli et al., 2006).In turn, this allows us to interpret the stadial/interstadial changes in sodium concentrations in the ice cores as qualitative indicators of the extent of the sea-ice cover in the North Atlantic.
In combination, the records of water isotopic composition, annual layer thickness and Na + and Ca 2+ concentrations in the ice allow us to study the phasing between changes in local temperature and precipitation on the Greenland Ice Sheet, North Atlantic sea-ice cover and dust deflation from the Central Asian deserts, respectively.To quantify this phase relationship we employ a probabilistic model of the transitions to determine their individual start and end points as well as the uncertainties of these points.The model describes the stadial-to-interstadial transition as a linear (in case of 18 O) or exponential ( :::::::: transition ::: (i.e.: a ::::: linear :::::::: transition ::::: fitted :: to ::: the ::::::::::::: log-transformed :::: data for all other records) transition between two constant levels, accounting for the intrinsic multi-annual, auto-correlated variability of the proxy records using an AR(1) noise process.This description of the transition is similar to the one used by Steffensen et al. (2008), using ::: who ::::: used a bootstrapping parameter grid-search algorithm to determine the ramp parameters and their uncertainties (Mudelsee, 2000).Instead ::::::: However, here we use probabilistic inference to determine the parameters of the ramp as well as those describing the multi-annual variability of the records and their uncertainties, conditioning on the data.Inference is performed on the model using an ensemble Markov Chain Monte Carlo sampler (Goodman and Weare, 2010;Foreman-Mackey et al., 2013) to obtain posterior samples for all parameters of interest as well as the parameters of the AR(1) noise model.In this way : , any correlations between parameter estimates and their uncertainties are transparently accounted for.More details on the inference and the mathematical description of the probabilistic model can be found in the Appendix.In addition to the start-and endpoints ::: end :::::: points we estimate the temporal midpoints ::: mid :::::: points of the transitions, as they are less influenced by the multi-annual variability that dominates the uncertainties of the timing estimates.However, here we focus on the interpretation :: of ::: the phase relationship of the parameters at the onset of the stadial-to-interstadial transitions to constrain the causal relationships of the trigger of the DO-events.
The transition model is individually applied to approximately 500 yr sections of the data at the highest constantly available temporal resolution (i.e.number of years per observation) around the stadial to interstadial transitions in each of the records on their respective time scales (Svensson et al., 2006(Svensson et al., , 2008;;Rasmussen et al., 2014).The exact width of the window was adjusted to account for gaps in the data and the onset of other transitions within the 500 yr window.Because of the thinning of layers due to ice flow, the time resolution decreases from 2 yr for the onset of the Holocene to 3 yr at the onset of GI-17.2 for the NGRIP CFA data, from 2 to 4 yr for the NEEM CFA data and from 4 yr at the onset of the Holocene to 7 yr for the NGRIP isotope data.The annual layer thickness data was down-sampled to match the NGRIP CFA data resolution.Note, that using the overall lowest available temporal resolution (7 yr) for the analysis leads to practically the same results, albeit with slightly larger uncertainties.It is also worth noting that the largest source for uncertainties in our estimates stems from the multi-annual variability of the proxy records that cannot be alleviated by higher resolution records.
An example of one of the fitted transitions, the onset of GI-8c in Ca 2+ , is shown in Figure 2 alongside marginal posterior distributions for the onset, midpoint ::: mid ::::: point and end of the transition.The absolute timing estimates for each proxy are used to infer the leads and lags between the different records, propagating all uncertainties.In the following, all estimates from the probabilistic inference are given as marginal posterior medians and their uncertainties as marginal 90% credible intervals.
As pointed out above, the cross-dating uncertainty between the two ice cores does not allow for sufficiently precise inter core-comparison :::::::: inter-core :::::::::: comparison of the absolute timing of the transitions, however the comparison of the lags between aerosol records of on ::: one ice core to those of another are not effected ::::::: affected by these uncertainties.Under the assumption that the DO-events show the same imprint in both ice-core records, : this enables crosschecking of the results between :: the : two records.
Note that this implicitly assumes that the timing differences for all insterstadial ::::::::: interstadial : onsets in the records ::::::::: parameters investigated here are the result of the same underlying process, : or in other words : , are similar between the interstadial onsets.This assumption differs from the assumption used in other studies of relative phasing of southern and northern hemisphere climate and changes in precipitation source regions during the last glacial (? Markle et al., 2016;Buizert et al., 2018) :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: (WAIS Divide Project Members, 2015; Markle et al., 2016;Buizert et al., 2018).In these studies, the stacking of the climate events assumes that the complete progression of the climate event is a realization of the same underlying process which is a wider and more restrictive assumption as it encompasses the whole transition and not only its onset.The combined estimates are calculated for the data from the two cores separately so that the consistency of the timing between Na + and Ca 2+ between records can be tested.Probability density estimates for the lead of the other parameters relative to Na + at the transition onset, midpoint and endpoint ::: mid :::: point :::: and ::: end ::::: point are shown in Figure 4.They clearly show thaton average , ::: on ::::::: average, both the reduction in terrestrial aerosol concentration as well as the increase in annual layer thickness precede the reduction in sea-salt aerosol for all stages of the transition, whereas no significant lead or lag is identified between 18 O and Na + .In the combined estimate, Ca 2+ concentrations start to decrease 7 +6 6 yr before Na + in the NGRIP record and 8 +5 5 yr in the NEEM record, where the error margins refer to the 5 th and 95 th marginal posterior percentiles.They reach the midpoint ::: mid :::: point : of the transition 9 +3 2 yr and 8 +4 3 yr before Na + and the end point 6 +6 5 yr and 4 +6 5 yr, respectively.

Discussion
Using much lower resolution mineral dust and ice 18 O data from NGRIP, Ruth et al. (2007) reported coinciding onsets within 5 to 10 years and a combined lag of 1 ± 8 yr for GI-1 through GI-24 demonstrating the close connection between the Asian and North Atlantic climate during all of the DO events.In light of the low temporal resolution of the data used (0.55 m corresponding to 7-98 yr or 38 yr on average) and the small timing differences that emerge in the combination of the individual GI :::::::: interstadial : onsets presented here, these results ::: the ::::: results ::: of ::::::::::::::: Ruth et al. (2007) are compatible with the outcomes presented here.Our results are also in good agreement with the more detailed studies of the onset of the Holocene , ::: and GI-1 (Steffensen et al., 2008) and GI-8 ::::: GI-8c : (Thomas et al., 2009).Similar to the presented study, Steffensen et al. (2008) andThomas et al. (2009) , (using high-resolution data from the NGRIP ice core, ) : also infer slight leads of changes in terrestrial aerosol concentrations and accumulation rate and moisture sources ahead of the changes in marine aerosols and local temperature in good agreement with our results from the respective DO-onsets.These results are also overall in good agreement with the other transitions investigated here.As the study presented here now covers 19 warming events starting from 60 ka ago and independently investigates the phasing in two ice cores, this adds significant evidence to the initially inferred phasing.
Two studies into Greenland (Schüpbach et al., 2018) and Antarctic (Markle et al., 2018) aerosol records have recently shown that a large part of the variability on centennial and millennial time scales in these records can be explained by the tight coupling between the hydrological cycle and the aerosol transport to the ice sheets.This coupling is a direct result of the efficient removal of aerosols during transport by precipitation scavenging leading to a negative correlation between en-route precipitation and aerosol concentration in the ice.On the short time scales investigated here, this tight coupling breaks down.This has previously been indicated for Antarctica by the lack of coherence at shorter periods (Markle et al., 2018) and here by the lack of synchronous changes between snowfall, water isotopic composition and/or aerosol concentrations during the onset of the warming events.This can be explained by the observation that under glacial conditions, models show that the interannual variability of the accumulation is governed by changes in the frequency of the advection of marine air masses to the Greenland ice sheet ::::::::::::::::::::::::::::::: (Pausata et al., 2009;Merz et al., 2013).Furthermore, the overall reduced amount of snow fall :::::: snowfall : during the glacial is a result of generally decreased moisture advection (Merz et al., 2013).For the stadial/interstadial transitions, these model results and the lack of covariance :::::::::: synchronous ::::::: changes : in the data presented here suggest a large influence of atmospheric dynamics on the change in snow accumulation as well, : .:::: This :: is further supported by other proxy evidence (Kapsner et al., 1995) :::: from :::::: central ::::::::: Greenland :::::::::::::::::: (Kapsner et al., 1995) ::: and :: is :::::::: similarly :::::::: observed :: in ::::: West ::::::::: Antarctica :::::::::::::::: (Fudge et al., 2016).An increase of the advection of air masses from lower latitudes would result in an increase of the snow accumulation as more moisture is transported to the ice sheet but in unchanged sea-salt concentrations in the precipitation and subsequently the ice core, as these air masses still carry the same amount of sea-salt aerosol.This leads to a decoupling of changes in the sea-salt aerosol deposition flux vs. : the concentrations in the ice i.e., the deposition flux per event stays constant but the number of :::::::::: precipitation : events increases.In turn, this means that the covariance of snow fall ::::::: snowfall : on the sea-salt deposition that is apparent on centennial and millennial time scales (Schüpbach et al., 2018;Markle et al., 2018) is less important at the interannual resolution :::: scale, where atmospheric dynamics dominate over thermodynamic changes.This view of changed precipitation frequency at the very beginning of the DO-onset is supported by the apparent asynchrony of snow fall ::::::: snowfall : and 18 O.Similar to Na + ,
Mineral dust, like sea-salt aerosol, is deposited very efficiently during snow fall ::::::: snowfall events because small amounts of precipitation remove most of the aerosol load from the air column (Davidson et al., 1996;Zhang et al., 2014).That means that an increase in snowfall frequency would not effect :::: affect : the observed dust concentration in the ice, despite the changing snow accumulation.As pointed out by Schüpbach et al. (2018), on centennial to millennial time scales, changes in the Central Asian dust source strengths :: by : a ::::: factor :: of ::::::::::::: approximately : 4 : are needed to explain the complete amplitude between stadial and interstadials observed in the ice cores, as further supported by downwind sediment records (Porter and Zhisheng, 1995;Jacobel et al., 2017).The deflation of the dust from the Central Asian source regions is strongly dependent on location of the westerly jet (Nagashima et al., 2011;Roe, 2009) that under current conditions co-varies with the East Asian Summer Monsoon on seasonal time scales (Schiemann et al., 2009).During the glacial conditions : , the westerly jet was located further south than today leading to more frequent or even permanent conditions for dust emission (Chiang et al., 2015;Nagashima et al., 2011).As low-latitude proxy records indicate northward movements of the ITCZ during the DO events (Wang et al., 2004;Deplazes et al., 2013;Wang et al., 2008;Yancheva et al., 2007), synchronous within approximately 180 yr (Adolphi et al., 2018), we hypothesize that the accompanying change in atmospheric circulation causes a reduction of the dust deflation during interstadials.Because the reduction in mineral dust aerosol concentrations in the ice core occurs before the reduction in sea-salt aerosol, the data implies that this change in atmospheric circulation happens before the reduction in the sea-ice cover in the North Atlantic.This is further supported by the observations by Steffensen et al. (2008), that show ::: fact :::: that changes in deuterium excess at the onset of GI-1 and GI-8 :::::: GI-8c, :::::::: indicating ::::::: changes :: in ::: in ::::::: moisture :::::::: transport :::: from ::: the :::::: North ::::::: Atlantic :: to ::::::::: Greenland, : occur after the onset of the transition in mineral dust aerosol :::::::::::::::::::::::::::::::::::: (Steffensen et al., 2008;Thomas et al., 2009).Antarctic records of deuterium excess that indicate changes of moisture sources also show a rapid shift at the onset of the Greenland interstadials interpreted as a change in the Southern Hemisphere Westerlies :::::::: westerlies : explaining part of the signal observed in multiple Antarctic
Even though the inferred timing differences carry large uncertainties for single events, they are overall consistent between the two investigated ice-core records and over all DO events indicating their robustness.The qualitative agreement between the DO events allows for the combination of the timing differences from the individual DO events to estimate a better constrained succession of the environmental changes during the onset of the DO events.Both the initial increase in snow accumulation and the reduction of Ca 2+ occurs about a decade before the reduction in Na + and the increase in local temperature.Taken at face value, this sequence of events suggests that the collapse of the North Atlantic sea-ice cover is not ::: may ::: not ::: be the initial trigger for the DO events and indicates that synoptic and hemispheric atmospheric circulation changes started before the reduction of the high-latitude sea-ice cover that ultimately coincided with the Greenland warming.This ::: The progression of environmental changes revealed in the Greenland aerosol records provides an essential benchmark a ::::: good ::::: target : for climate models that aim to :::::::: explicitly :::::::: modeling :::: both ::::: water :::::: isotope ::: and ::::::: aerosol ::::::: transport :::: that ::: aim :: at ::::::::: transiently : simulate DO events.However, using ::::: Using our ice-core data, we can only derive the relative timing of changes in the atmospheric circulation and the North Atlantic surface conditions related to the DO-warming.Our results cannot provide a direct constraint on the increase in AMOC which may potentially precede changes in sea-surface conditions and atmospheric circulation.Recent studies using marine sediment records (i.e. Henry et al., 2016) ::::::::::::::::::: (e.g. Henry et al., 2016) suggest a centuries-long lead of AMOC changes, however the response time of the proxies, their limited resolution and dating precision may not allow an unambiguous answer yet.

P3 L27 :
All data ARE shown. . .(the word "data" is plural not singular) Updated accordingly.

Figure 1+3 :
Figure 1+3: Please consider plotting the age axis in reverse (so time goes from left to right).This is what much of the paleoclimate literature is moving towards.That is also what Figure 2 uses.
Bigler, M. et al. (2006), Aerosol deposited in East Antarctica over the last glacial cycle: Detailed apportionment of continental and sea-salt contributions, Journal of GeophysicalResearch, 11110.1029/2005JD006469.
lines 11: I think this should be "...a coinciding..."Not "...an coinciding..." Corrected, thanks.Page 11, lines 21-22: How big are the additional changes in source need compared to the total variability across a DO event?

Figure 1 .
Figure 1.Investigated records from the NEEM (a, b) and NGRIP ice cores (c-f).(a-d : a, : d) show the new aerosol records of Ca 2+ (a, c) and Na + (b, d) alongside the NGRIP layer thickness ( , : :: (e) and ice 18O records (f).All records are shown in :: as decadal averages on their respective version of the GICC05 age scale relative to 1950(Svensson et al., 2008;Rasmussen et al., 2013).The vertical lines mark the

Figure 2 .
Figure 2. Example of a fitted ramp for the Ca 2+ data from the onset of GI-8c.The top panel shows the data together with the marginal posterior median of the fitted ramp (thick black line).Shaded areas indicate the marginal posterior 5 th and 95 th percentiles.The time is given relative to the nominal transition as given in Rasmussen et al. (2014) with negative values meaning before the transition and positive after.The bottom panel shows the marginal posterior densities for the onset, midpoint ::: mid :::: point and endpoint :: end ::::: point of the transition (from left to right : to ::: left).

Figure 3 .
Figure 3. Timing differences for the individual interstadial onsets.: : (a) shows the NGRIP 18 O record, (b-d) timing difference of NGRIP Ca 2+ , layer thickness ( ) and 18 O records as well as the NEEM Ca 2+ record relative to the transition in the respective Na + records at the onset (b), midpoint ::: mid :::: point : (c) and end (d) of the transition.No timing results are given for transitions where there are data gaps in one of the necessary data sets.Error bars show the marginal posterior 5 th and 95 th percentiles and the symbol the marginal posterior median.Note the different axis scaling for the start, middle and end points.All inferred absolute timings, ::: the ::::::: transition :::::::: durations and ::: the lags relative to Na + can be found in the Supplementary Material.

Figure 4 .
Figure 4. Combined probability density estimates of the lag of the Ca 2+ , layer thickness and 18O transitions relative to the respective point in the Na + records for the onset (a), midpoint :: mid :::: point : (b) and end (c) of all transitions from stadial to interstadial for the two cores.