The paper titled “Magnitude, frequency and climate forcing of global volcanism during the last glacial period as seen in Greenland and Antarctic ice cores (60-9 ka)” is the first effort to reconstruct the volcanic history of the last glacial period, utilizing six bipolar ice core records. Results from this work will provide an essential forcing scenario for model simulations and climate investigation of this period. I recommend publication of this work, after careful consideration of the following questions and suggestions:

1. Section 2.2, please explain briefly the choice of different filter length, i.e., 45yr for the Antarctic cores and NGRIP while 181yr for GISP2 & NEEM.
2. Also in section 2.2, it is a bit confusing about which parameter was used to measure the volcanic signals, the depth or the yr. Based on the text, the background was filtered by window length indicated by year, while the duration of the volcanic signal was indicated by depth.
3. Section 2.4 please provide some description of how the manual correction was performed. For example, what are the resolutions for the ECM and DEP records? How many or what percentage of the signals have been corrected. It would also be great to give an example of how the correction was done.
4. “A sulfate deposition of 20 kg km^-2 corresponds to half the Greenland deposition from the 1815 AD Tambora eruption, thus refers to quite large events in terms of total sulfur injections into the atmosphere.” Please explain briefly how was the 20 kg km^-2 (and also the 10 kg km^-2 for Antarctic) cutoff line estimated. And was the 40 kg km^-2 1815 AD Tambora deposition corresponding to the average deposition from the three Greenland cores?
5. Section 2.6 Please provide more details on the SVM model. For example, what are the requirements, the pro. and cons of the model in this particular application. What validation had been done on the model performance? For example, taken one of the 21 eruption signals used to train the model out from the analysis, could the model accurately simulate its location?
6. Ln 305-307, the comparison between IC and CFA records in Fig 2e needs further demonstration. For example, what is the exact meaning of “very large uncertainties”? What is the implication of the uncertainties on the interpretation of the “face value”?
7. Section 3.3 Please explain why borrow the bipolar eruptions from the previous studies, rather than estimate a list using the results from this study.
8. Ln 532-534 In my understanding, the VEI list is a discrete (i.e., it is not a complete but continuously evolving) reconstruction of historical volcanism based on geological investigation. So I am not sure it is appropriate to directly compare the event frequency from the ice-core-based reconstruction (which is assumed to be continuous) with that from geological investigation.
9. Section 4.4 Please explain why do some events have forcing estimation in a range (for example, the #2 largest signal has forcing ranging from 17.8 to 176.5 W/m2), while others have finite forcing estimation (for example, the #3 largest signal has a forcing estimation of 82.8 W/m2)? If it was due to the number of ice cores available for signal extraction, this should be clarified.
10. Is there any reason why the authors choose Tambora & Samalas for the magnitude comparison of #3 and #4 events, respectively?
11. Please add captions for the supplementary figures.

This paper represents a huge amount of work to identify the deposition of volcanic sulfate in the polar regions and to try and draw conclusions about the occurrence of volcanism over the last 60 kyr. I really applaud the effort, and the attempt to draw large-scale conclusions from it. I have a lot of relatively minor comments on the paper (which put it somewhere between minor and major revision), and I think the authors could have made some better choices in the way they treated the data. However I accept that they made mainly reasonable choices and so I do not propose to insist on any significant reanalysis of the data – just in places an extra sentence is needed to discuss the choices made.  I found some of the messages that end up in the abstract too strong given the nature of the data and I will comment on those in the text, expecting the authors also to address them in the abstract.

A final point is that the paper is extremely hard to follow – while I appreciate the need for a lot of supplementary material, the fact that the figures in the main text and the supplement seem to be called in almost random order is very unhelpful, and I suggest the authors renumber the figures to provide a more logical flow.

Before I give a detailed set of comments, I should add that (as I think several of the authors are aware) I have been involved in a similar exercise, using only the EDC ice core, but for a period of 200 kyr. This has been presented a few times and has been (in a paper involving also terrestrial and marine data) under review for a considerable time. My comments should therefore be taken with the knowledge that I have used slightly different methods and criteria but am not seeking to impose the same methodology on the authors here.

Detailed comments:

The paper will need a thorough copy edit as the English at times is awkward and occasionally hard to follow.

Line 165. I think the authors have misunderstood the MSA correction. This had to be done for WAIS Divide because ICPMS measures elemental S, and therefore the values obtained are (sulfate + MSA). This is not the case for any of the FIC/IC methods, where sulfate itself is measured, so no correction would be required and no correction should have been made for any other sites. Please clarify this in your text.

Line 188. I found myself a little confused about the volcanic detection threshold. I understand the idea behind the use of RRM and RMAD. But then in the end you just use a threshold of 10 or 20 kg km-2 so I am not sure what purpose the statistical threshold serves.

Equation 1: please be careful to define the units when describing this equation.  I assume concentration is in ppb or ng/g (as in the figures), and D is in m ice equivalent (this is important as if you had included the top 100 m where density is less than that of ice, then the equation is not correct if in real depth), and 0.917 is in g/cm3. By fluke, after cancelling factors of 10 to change the units this does indeed work out as kg km-2, but the reader needs that to be stated.

Line 209. Please also be careful to define what T is, ie layer thickness/original layer thickness. It’s confusing otherwise because you talk next about a 60% reduction of the layer thickness which is equivalent to a T of 0.4. Perhaps that is the best way to deal with it, “60% reduction of the layer thickness (T=0.4)”.

Supp table S1: Please explain what the column “Name” is: I assume it’s the age model used.  But I’m surprised because the standard for a 60ka period is GICC05. I appreciate you may have to use a model to derive a smooth acc rate, but then you need to explain the relationship between ss09seabm1 and GICC05 as few of your readers will have heard of the former.

Fig S1 and caption. Please clarify that in part c the y-axis “thinning function” is equivalent to “T” in the text as I don’t think you ever explain that in the text. Also in the caption, what does “thinning file” mean? The second sentence of the caption doesn’t make sense, please re-word.

Line 220. This is wrong. The volcanic sulfate flux calculation does not assume anything about wet or dry deposition: it is simply the product of concentration and snow acc rate. It’s true that the flux is affected by whether there is significant wet deposition in addition to the dry deposition and this is a point worth making, but the flux as calculated is definitely correct.

Table S4, column O, kg not Kg please.

Table S5. Why is N/A written for WD2014 ages below 30 ka. WD2014 extends from the surface to 68ka even if it changes from layer counted to methane-tied at 31 ka.

Tables S4 and S5: Are the values given for individual sites before or after the rescaling of Antarctic concentrations. I think before – if not then it’s hard to understand the averages you give for the Antarctic. Please clarify this in the captions, ie that the individual values are as measured while the average is after rescaling.

Line 241. What matters is of course not the value of the sulfate background but its variability. I think this is what you mean but from this wording it isn’t quite clear. I think you need to explain your 10 and 20 kg m-2 criteria better – how is it derived? In my own work we estimated the “negative” peaks to understand how big a deviation from the median could be generated in a given section from noise alone, and then we used that to establish a threshold (which was 20 kg km-2 for the 200 ka record). So I think your thresholds are probably fine, but could be explained more clearly.

Line 255. The Antarctic plots shown in Fig S2 are really not impressive with very small r^2. They don’t seem at all a good basis for the rescaling you do. While I accept there may be higher fluxes for a given volcano at WD because it receives more wet deposition (this is the only reason the accumulation rate is relevant), the values you derive here have a huge uncertainty which you have not apparently propagated into your 26% error estimate (line 265). We know there is a huge uncertainty anyway because of local variability in deposition (as shown eg by Gautier et al (2016). I suggest you do an alternative error estimate where you use the rescaled data from the 123 Antarctic eruptions where you have 3 sites and find the average 1-sigma and/or 2-sigma between the 3 sites. This would give you an alternative way to estimate the uncertainty in the values you end up with.  26% is certainly way too small an error in the cases where you have only 1 site.

Line 268. I am familiar with the paper by Marshall (not Martshall) but I don’t know what you mean by “The volcanic sulfate deposition in Greenland and Antarctica shows a distribution pattern (Table S5 and Table S6), similar to that derived from the aerosol-climate modeling of volcanoes over past 2500 years”. I can believe that you are assuming they show such a distribution and that this is how you decide on the latitude, but you have no evidence to say that they do show such a pattern.

Section 2.6 and Fig S3. I apologise for my lack of technical knowledge but this section is not comprehensible to someone coming at it new. I think you need to explain what you are trying to achieve here (which I think is to use the Greenland/Antarctic ratio of known eruptions to classify the latitude of unknown bipolar eruptions). Then please try to explain SVM better. To anyone not in the know Fig S3 b and c cannot be understood.

Line 282, I think you mean Fig S10, not S9.

General comment: this reminded me that the calling of figures, especially supplementary ones, in apparently random order is really confusing. Please sort this out.

Line 360. I don’t really understand your statement that “the layer thinning becomes stronger which makes it increasingly difficult to detect smaller eruptions signals”. If that is really so then your threshold of 10 or 20 kg km-2 is too low. The whole point of it should be to ensure that detection is the same throughout.

Section 3.4. The previous comment raises a more general comment on this section. It’s really obvious that you shouldn’t use cores where the resolution is worse than the expected peak width. However it’s difficult to know what that means: in the case of EDC, if thinning were the only thing happening then it would be hopeless to detect eruptions reliably in the early parts of the record with resolution 5 years. However it turns out that the diffusion at Dome C keeps the peak width quite constant with age, and therefore allows resolution of peaks where the age resolution is 5 years (because the peaks at this depth have diffused to 10 or 20 years wide). However I’d be very surprised if you can expect to resolve most volcanoes, even large ones, in GISP2 with the stated resolution in Table S1. For NEEM, given the variable depth resolution, I suspect the age resolution is normally much better than 10 years (I suggest checking this again), but if it were 10 years again I think detection would be hopeless. I think you need a somewhat deeper discussion here, and also to consider whether it is worthwhile to contaminate a great dataset by including data from sites where detection must be severely degraded.

Fig S5 and line 375. Please be careful here. I believe the Rasmussen analysis makes some sense for the kind of CFA set up used by Bigler where there is a large reaction cell that is being continually mixed leading to quite a large mixing volume. In such a setup there can be a nominal data resolution of mm which is not real. The FIC setups are a bit different as they inject a slug of liquid every 2 cm (or whatever each system uses) so the resolution could never be better than that in any sense. I’m not sure your analysis really makes that distinction.

Line 389. You call Fig 4 before Figs 2 and 3.

Line 419. I found this very confusing. I am looking at Fig 2 where you show the Holocene, deglacial, stadial and interstadial. You then state that the number of eruptions is higher in cold periods when the figure shows the most in the deglacial and really no difference between the glacial and Holocene. Eventually I realised that you are only talking about stadials versus interstadials within the glacial. So firstly please make this clear throughout this section. But really is that statement true at a significant level: perhaps so in the 1 year data in Fig S6, but the reader isn’t looking at that, they are seeing Fig 2. You need to be much clearer that you are using Fig S6 to discuss whether the difference is an artefact and that the conclusion, as illustrated in Fig 2 is that it probably is. I realise you do say this in the end, but by pointing at Fig 2 initially you leave the reader who doesn’t take the trouble to look at the supplement lost.

Line 488 (2-0) Do you mean 2-0 ka? When I look at Fig 2, I don’t really see a 34% increase. At least for the small (blue) category, there is no change. Perhaps describe this more precisely as to what you mean.

Line 535-544. I don’t find this comparison of the ice core eruption rates against Lameve data helpful. Firstly you should also look at the analysis by Rougier et al (Rougier, J., R. S. J. Sparks, K. V. Cashman, and S. K. Brown (2018), The global magnitude–frequency relationship for large explosive volcanic eruptions, Earth planet. Sci. Lett., 482, 621-629, doi:https://doi.org/10.1016/j.epsl.2017.11.015.). It is very clear in that paper that the community is well aware that Lameve under-reports so I don’t really see the value of the comparison you are making.

I am unconvinced by section 4.4 Given the huge uncertainties, especially for eruptions where there is only 1 core represented in one of the ice sheets, the values are so uncertain that trying to call out individual positions in the medal table seems a bit pointless. At the least please add a column where you state how many cores are represented in each ice sheet. But  I feel this section of the paper is given too much weight and should be shortened. Remember also that you certainly show no eruptions that are bipolar between 16.4 and 24.5 ka, but this is not because they were absent but only because the tie points are not good enough to know whether they are bipolar. Given this kind of issue I think the league table, while it can be included should be downplayed and put into context.