Comments on: Using reduced-complexity volcanic aerosol and climate models to produce large ensemble simulations of Holocene temperature

Verkerk and coauthors combine reduced-complexity volcanic aerosol (EVA_H) and climate (FaIR) models to simulate the global mean surface temperature (GMST) response to volcanic eruptions over the last 9,000 years (6755 BCE to 1900 CE).

To assess the robustness of their simulations, the authors compare their estimates for the 14 largest eruptions between 1250 CE and 1900 CE with numerous climate reconstructions (Schneider et al., 2015; Wilson et al., 2016; Guillet et al., 2017; Pages2k, 2019; King et al., 2021). The discrepancies between the new simulations and climate reconstructions are notably smaller than in previous studies.

The authors address an important topic. The paper is well-written, well-structured, and easy to follow. The figures are clear and informative. And the authors have made all their simulations publicly available.

The methodology section summarizes well the approach taken by the authors, including the forcing datasets used for the new simulations, the paleo-reconstructions and the climate simulations employed to compare the new results.

Additionally, they acknowledge the limitations of their approach, particularly the Holocene temperature conundrum, which is also apparent in their ensemble simulations of Holocene temperatures.

The authors emphasize the need for future products based on reduced-complexity models to include seasonal and regional outputs, which would be highly valuable for the paleo community.

I appreciated reading the manuscript and, overall, have very few comments to offer. I recommend the paper for acceptance, as I think the new product provided by the authors represents a valuable resource for the paleo community studying past volcanic eruptions. However, I do have one minor suggestion for the authors to consider.

Main text:

• Comparing simulations with instrumental data: Pushing the simulations beyond 1900 CE would have been a great addition. Extending the simulations into the 20^th century would allow direct comparisons with instrumental data for eruptions such as the 1902 (Santa María), 1912 (Katmai/Novarupta), 1963 (Agung), and 1991 (Pinatubo) events. They could help validate the accuracy of the simulations.

Have the authors considered the possibility of comparing the accuracy of their simulations not only against climate/data assimilation reconstructions but also against instrumental datasets, such as the Berkeley Earth Surface Temperature (BEST) dataset? The BEST dataset offers two products that might be of interest: one estimating GMST since 1850 and another providing annual temperature estimates since 1750 (land-only).

Using these datasets could allow the authors to compare their simulations for the 1815 Tambora, 1831 Zavaritskii (Hutchison et al., 2024), and 1883 Krakatau events with “real” temperature observations. Additionally, the Laki eruption might also be investigated, assuming the instrumental records used by BEST are sufficiently dense to represent a reliable global average (which I am not entirely certain about).

• Line 130: Change Hutchison et al., in review to Hutchison et al., 2024

Supplementary Material

• Line 60: “Table S3: Integrated response of the superposed epoch analysis (Error! Reference source not found.d).” There appears to be a reference issue here that should be corrected.

The authors use a reduced-complexity climate model to obtain an ensemble of simulations of the last 9000 years with different volcanic forcing for the purpose of quantifying volcanic forcing uncertainty, and to put them into the context of proxy reconstructions. This is a nice and well written paper, but I have some concerns about their conclusions. I hope the authors will accept my constructive criticism which is meant to increase the novelty of the manuscript and to strengthen their results.

Major remarks:

(1) My main point of concern is the fact that the NH climate has been simply scaled to match the observed record for the 20^th century. While the authors have noted in the discussion that this approach is simplistic, two of their three major conclusions rely heavily on their NH temperature simulation and thus on a simple scaling approach.

The problems I see are:

• Scaling suppresses the amplitude of the record and thus the response to volcanic eruptions.
• However we do not know if the same scaling applies during periods of volcanic activity.
• In fact volcanic activity is almost negligible within the 167 yr record considered for the linear regression.
• Thus the scaling approach is not validated for representing volcanic years, however this is the main purpose of the paper/

In short, we may see a better fit between the here presented simulations and proxy reconstructions due to the suppressed amplitude, but this can likely be an artefact of the scaling approach. Their results therefore need to be strongly caveated, which significantly reduces the novelty of this manuscript. I would suggest to either obtain a better way of quantifying NH climate, or to put a larger emphasis on their global results.

In order to provide a better context for their results, they could also repeat their model simulations for volcanic forcing using EVA instead of EVA_H and compare the results. This would provide a better context for comparisons with Luecke et al. 2023 and with reconstructions. If they can thus prove that the reason for the better performance results from the use of EVA_H compared to EVA, and not from their scaling, this would justify their conclusions. If the different results are a result of the superiority of the FAIR model over the model used by Luecke et al., the authors need to find a better way of estimating NH climate.

(2) Another point of concern is that dating uncertainty has not been taken into account. For the purpose of comparison with proxy reconstructions, this is a very important source of uncertainty, and in particular plays a major role for presenting superposed epoch analyses. Since the SEA is one of the major results presented, I would suggest to repeat the SEA but to perturb the sampled eruption year and thus get a measure of how dating uncertainty would change the amplitude of their simulations.

(3) I am also intrigued that the authors have provided a full 9000 year simulation, however concentrate largely on their discussion of last Millennium climate. Would it be possible to extend their discussion of Holocene climate? As I understand this approach of using a reduced complexity model and a large forcing ensemble is completely novel for Holocene climate. I would like to have a clearer idea what we learned from this experiment!

(4) I am not sure how internal variability is accounted for. It is mentioned a few times but nowhere explained in detail (or has slipped my attention).

(5) Lastly, I understand the authors have used different implementations of land-use forcing. As far as I can tell no results have been shown or are being discussed anywhere. Can you please clarify what the purpose of this experiment was?

Further comments:

• 104 ff. How is the eVolv2k record different to the ensemble described in Luecke et al. 2023? What is the added benefit for creating a new ensemble from scratch? I think this needs to be more clearly brought out in the introduction. Is the novelty the length? Or the new approach?
• 3: M has not been introduced. I understand that the reference has not been published yet. But this needs further justification or context (i.e. where do the scaling factor and exponent come from?).
• Why is timing uncertainty not included???? This is a key uncertainty when comparing to proxy data, and in particular heavily affects the SEA??
• L 265: ‘EVA_H’s empirical nature’ – In what sense is this empirical? Is this the same for EVA? Please clarify.
• L 267-269: So caution is needed for interpreting both eruptions larger than Pinatubo and much smaller than Pinatubo. How would this bias the results?
• L 268-269: Please can you add a quantifying statement to this.
• L 290: anthropogenic aerosol and ozone forcings – this is new to me as has not been mentioned anywhere before? Should be added to the section about forcing. Also please add reference.
• L 320: I do not understand how the design accounts for internal climate variability!
• 370 ff. Given the extensive discussion of the MCA-LIA difference it would be worth adding a boxplot figure showing the differences between the reconstructions.
• Figure 6: I wonder what cause the deviation between your simulations and Luecke et al. Is this simply a result of the anomalies and a mismatch in the 1850-1900 period or is this the result of different forcing or model differences?
• Section 5.: The fact that the MCA-LIA difference is much smaller in the simulation compared to the proxy reconstructions could also be an artefact of tree-ring spectral biases, in particular RW type proxies (which also explains why Schneider does not show this difference). Please add this to the discussion and reference Luecke et al. 2019
• Section 5.4 has a large overlap with the conclusions and could be merged.

Minor remarks:

• L18-19: What about dating uncertainty?
• L 22: So no anthropogenic aerosols?
• L 25: is this averaged over 9000 years?
• L 32: “a relatively cool period in climate reanalyses” – Wording is unclear. Do you mean the cool period is found in climate reanalyses or is the cooling not found in reanalyses?
• L 46: tree ring data is the only proxy with reliable annual resolution for the last Millennium.
• L 47: Cite Luecke et al. 2021 for seasonal bias
• L 53: Cite Luecke et al. 2023 for volcanic forcing uncertainty
• L 72: Do you mean “…we refer to ‘reduced-complexity models’ as idealised models”?
• L 126: I’m not sure if this is due to the draft setting, but the formatting here is really off (also in the following equations)
• L 127: asym is in math mode (italic) but should be text mode (upright) in eq. (1)
• L 131: Is it worth quantifying the uncertainty for identified eruptions? I assume there is significant uncertainty associated with those estimates of plume height. It would be interesting to know how this compares with unidentified eruptions uncertainty.
• L 131: Is this correction about identified eruptions? This is a bit hard to follow for anyone unfamiliar with isopleth maps.
• L140/eq. 3: What is M? This has not been introduced.
• L 141: confusing wording: Fig. 1f provides an overview but emission parameters are in SI? This is unclear. Also include reference to where in SI it’s shown.L 201-209: No remark but I really like this brief discussion of limitations of the proxy records.
• This discussion of the key differences to Luecke et al.’s study should have been discussed at least briefly in the introduction.
• L 258: ‘EVA_H accounts for …’
• L 274: put \alpha in math mode
• (4) put -\Delta gmSAOD into text mode
• L 275: It would be better to write this in text, i.e. ‘the relationship between gmSAOD and ERF’ otherwise it could be interpreted as gmSAOD minus ERF.
• L 276 Eq. (4)
• L 294-295: does this ensemble represent different values for climate sensitivity? Please clarify.
• L 327-328: Okay but uncertainty is also very large for these eruptions. So within the uncertainty range (assuming the gmSAOD uncertainty is symmetrical around the mean), more recent eruptions could in fact exceed them.
• L 333-334: Add commas to help flow: Volcanic injections with known eruption match, and thus better constrained latitude and altitude, have smaller forcing uncertainties…
• L 395: Wording is very clunky with the double negation, please clarify.
• L 487-490, L516-518: Strongly overemphasises the conclusions for NH climate here. Need better quantification/justification why the agreement is better, otherwise please caveat results.
• L 519-520: how can discrepancies be explained by internal variability? I really don’t see where this has been taken into account.

• Fig 1: I like the figure, especially fig. 1f which is a nice visualisation of eruption parameters. However I’d recommend reordering and start with 1f, which is (i) your most important figure and (ii) mentioned first in the text.
• 4b shows GMST for single forcing runs- does it show the ensemble mean for volcanic forcing? And just one implementation of land use. Can you clarify the choice?
• 5: it is striking that the model here show an ongoing warming trend and a lot more variability at multi-millennial scale than the proxy based data. What are the reasons for this? If Kaufman and Broadman 2023 suggest this could be from volcanic forcing and your results rebut this then it would be worth putting more emphasis on this.