Large uncertainty in volcanic aerosol radiative forcing derived from ice 1 cores 2

Large uncertainty in volcanic aerosol radiative forcing derived from ice 1 cores 2 Lauren Marshall, Anja Schmidt, Jill S. Johnson, Graham W. Mann, Lindsay Lee and 3 Ken S. Carslaw 4 5 School of Earth and Environment, University of Leeds, Leeds, UK 6 Department of Chemistry, University of Cambridge, Cambridge, UK 7 Department of Geography, University of Cambridge, Cambridge, UK 8 National Centre for Atmospheric Science, University of Leeds, Leeds, UK 9 10 Corresponding author: Lauren Marshall (lrm49@cam.ac.uk) 11 12 This is a non-peer reviewed preprint submitted to EarthArXiv 13


33
Explosive volcanic eruptions that inject large amounts of sulfur dioxide (SO2) into the stratosphere are 34 major drivers of natural climate variability on multi-annual to decadal timescales 1 . The SO2 is 35 converted to sulfate aerosol, which causes a radiative perturbation, or forcing, by scattering incoming 36 solar radiation, and leads to surface and tropospheric cooling. Reconstructions of volcanic aerosol 37 radiative forcing are therefore required to understand and attribute climate variability on millennial 38 timescales 2-4 and are used as input to climate model simulations 5,6 . Correct reconstructions of all 39 climate forcing agents are vital to understand and evaluate past temperature changes on global and 40 regional scales, to assess ocean heat uptake, which is critical for estimating future sea-level rise, and 41 ultimately to compare natural and anthropogenic drivers of climate variability. 42 size distribution 29 , the latitude and season of an eruption 30 , aerosol-cloud interactions 31-33 and model 73 configuration 1,28,34 . 74 Simulations of the last millennium as part of the Paleoclimate Modelling Intercomparison Project 75 phase 4 (PMIP4) will use the prescribed sAOD forcing timeseries 'EVA(2k)' derived by Toohey and 76 Sigl (2017) 18 , where the radiative forcing is calculated internally by each model. The spatial and 77 temporal evolution of the prescribed sAOD is based on a simple parameterized transport model, the 78 Easy Volcanic Aerosol (EVA) forcing generator 35 , which uses stratospheric SO2 emissions from the 79 eVolv2k reconstruction 18 . The EVA forcing generator is calibrated against the measured evolution of 80 stratospheric aerosol following the 1991 Mt. Pinatubo eruption and does not account for many 81 microphysical, chemical and dynamical processes that occur following a volcanic eruption. Because 82 many eruptions identified in ice-core sulfate records are unattributed, the eruption season and latitude 83 must be estimated or arbitrarily assigned, which introduces further uncertainty because these factors 84 affect the formation and transport of stratospheric sulfate aerosol and its deposition 17,22,36 . Eruptions 85 are assumed to be tropical if simultaneous sulfate signals occur in both Antarctica and Greenland 86 (bipolar deposition signals) and are assigned to January if the season is unknown. 87 The difficulty with any reconstruction of radiative forcing is that it does not scale directly with the 88 ice-sheet-deposited sulfate. The magnitude of the forcing (integrated over time) depends on the global 89 spread of the volcanic aerosol in the stratosphere, its lifetime, and the microphysical properties of the 90 aerosol (size, mass and number of the aerosol particles). All of these depend on the emission strength, 91 the altitude and latitude of emission and the eruption season, so-called 'eruption source parameters' 37 . 92 Consequently, for any observed ice core volcanic sulfate deposition there is potentially a very wide 93 range of 'eruption-realisations' (i.e. eruptions with different combinations of SO2 emission, eruption 94 latitude, emission altitude and eruption season) with a wide range of associated forcings. 95 Some attempts have been made to estimate the uncertainty in SO2 emissions derived from ice-core 96 sulfate composites by considering uncertainties in the ice-core composites themselves and in the 97 transfer functions 18 but the possible range in radiative forcing has not yet been properly quantified.  On average, the deposition of sulfate (SO4) on Greenland is higher than on Antarctica, with a 136 maximum of 148 kg km -2 deposited for an eruption at 79ºN occurring in January, with a SO2 emission 137 of 84 Tg. The maximum simulated Antarctica deposition of 65 kg km -2 occurs for a July eruption at 138 72ºS with a SO2 emission of 98 Tg. Lower deposition on Antarctica compared to Greenland was also 139 found in a previous study of tropical eruptions 22 , most likely due to stronger meridional transport in 140 the Northern Hemisphere (NH) and increased deposition because the NH is relatively more 141 dynamically active than the Southern Hemisphere (SH). In the SH the stronger polar vortex will 7 inhibit more of the poleward aerosol transport. Deposition on the ice sheets will also vary with SO2 143 emission magnitude given an increase in sedimentation as particles grow larger such that they may be 144 deposited before reaching the ice sheets, and stronger polar vortices arising from aerosol-induced 145 stratospheric heating 22,38 . 146 Deposition on Greenland is greater for tropical eruptions occurring in January (blue circles in Fig. 1a) 147 because more sulfate aerosol is transported to the NH via the Brewer Dobson Circulation (BDC), 148 which is stronger in the winter hemisphere. Similarly, both total SH deposition ( Supplementary Fig.  149 1) and deposition on Antarctica from tropical eruptions is greater if they occur in July (red circles in 150  Fig. 1), however ice-sheet 152 deposition varies between seasons, but is not consistently larger in either one. Differences in the ice-153 sheet deposition following eruptions in different seasons for mid-to-high latitude eruptions could be 154 dependent on the SO2 emission magnitude, injection height, and seasonal variations in stratosphere-155 troposphere exchange and sulfate aerosol deposition rates in the mid-latitude storm tracks 39 . Seasonal 156 differences may also arise due to internal variability. 157 There is also considerable scatter in the deposition values around zero for eruptions located at high 158 latitudes in the opposite hemisphere to the ice sheet, which has implications for the detection and 159 quantification of volcanic events. We find that the ice-sheet deposition time series are very noisy 160 because of internal variability, which can obscure or enhance the deposition ( Supplementary Fig. 2). 161 This is because the difference between the ice sheet sulfate deposition in the volcanically perturbed 162 simulations and in the control simulation (see Methods) is affected by volcanic sulfate deposition, as 163 well as by variations in the background tropospheric sulfate aerosol deposition caused by the effect of 164 the eruption on stratospheric and tropospheric dynamics (i.e. the control and perturbed runs 165 effectively behave like two meteorological ensemble members). Consequently, the amount of 166 background tropospheric-originating sulfate aerosol can be very different in each perturbed simulation 167 compared to the control climatology. Time-integrated deposition anomalies can even be negative 168 because the climatological deposition is higher than in the perturbed simulation, which represents just one possible realisation of reality (an example is shown in Supplementary Fig. 2). Our analysis of a 170 wide range of simulated eruptions highlights the inherent difficulty in detecting and quantifying 171 volcanic deposition anomalies in ice cores. The volcanic sulfate deposition signal on Greenland 172 becomes clear only for anomalies that exceed ~20 kg km -2 and for Antarctica when anomalies exceed 173 ~10 kg km -2 ( Fig. 1; grey lines). 174

Ice-sheet sulfate deposition predicted for thousands of eruptions 175
By replacing the UM-UKCA model with statistical emulators that describe how the ice-sheet 176 deposition varies with SO2 emission, eruption latitude and injection height, we can predict the sulfate 177 deposition for any eruption that we did not directly simulate (see Methods). The emulators, which are 178 built for January and July eruptions separately, enable us to examine the relationship between ice-179 sheet sulfate deposition and the eruption source parameters in unprecedented detail because the 180 emulated predictions for all possible eruptions in our three-dimensional parameter space describe how 181 deposition varies continuously ( Supplementary Fig. 3), effectively interpolating between the model 182 output of the simulations. We account for the influence of internal variability by adding a noise 183 variance term to the deposition anomaly during the construction of the emulators such that we do not 184 need to run conventional meteorological ensemble members for each eruption realisation (see 185

Methods). 186
The combinations of eruption source parameters that lead to deposition within a measured range can 187 be found by constraining the emulator predictions for all possible eruptions in our three-dimensional 188 space (see Methods). To illustrate the constraint procedure we take the emulated deposition and find 189 the source parameters of eruptions that lead to the ice-sheet sulfate deposition derived from ice cores 4 190 for the 1815 eruption of Mt. Tambora (Table 1) For an eruption occurring in January, the predicted deposition falls within both the Greenland and 208 Antarctica deposition uncertainty ranges only if the SO2 emission is greater than 73 Tg and the 209 latitude of the eruption is between 20°S and 49°S. The injection height remains unconstrained. For an 210 eruption in July, only eruptions with SO2 emissions greater than 81 Tg and with a latitude between 211 4°S and 59°S can produce deposition that matches both ice sheet constraints. To match the ice sheet 212 deposition for an eruption occurring at the latitude of Mt. Tambora (8°S), the eruption must occur in 213 July and the SO2 emission must be greater than 96 Tg. This emission is higher than the 60 Tg estimate 214 often used to simulate this eruption in climate models 40,41 , but closer to petrological estimates (73-91 215 Tg) 15 . We have built emulators only for eruptions occurring on 1 January and 1 July, whereas Mt. 216 Tambora erupted in April 1815. The predicted deposition is generally higher in Greenland for January 217 eruptions, and higher in Antarctica for July eruptions ( Supplementary Fig. 3) and deposition following 218 an April eruption would also differ given seasonally varying stratospheric transport of sulfate aerosol 219 and depositional processes 17 . Furthermore, previous simulations of the 1815 eruption of Mt. Tambora 220 using UM-UKCA with 60 Tg SO2 emitted at the equator showed that sulfate deposited on the ice 221 sheets was roughly half that derived from ice-core estimates 40 . This indicates that either the SO2 222 emission used in the model was too low or that a structural error exists within the model resulting in a 223 low bias in deposition. Background (non-volcanic) sulfate deposition simulated in UM-UKCA was 224 found to compare well to ice-core estimates. 225

Constraining volcanic radiative forcing for the ten largest bipolar deposition signals 226
We now examine the full range of eruptions that could lead to the ice-core-derived deposition signals 227 for the ten largest bipolar events in the last 2500 years 4 . Only two of these events have been 228 confidently attributed to known eruptions (1815 Mt. Tambora, which is the 6 th largest deposition 229 signal, and 1257 Samalas, which is the 2 nd largest deposition signal) 18 . For each eruption-realisation 230 retained in the constraint, we examine the cumulative sAOD and cumulative RF predicted by the 231 respective emulators (see Methods) to estimate the range in volcanic forcing that is consistent with 232 each deposition signal. 233 The 44 BCE and 266 CE eruptions deposit much more sulfate on Greenland than on Antarctica (a 309 ratio of 6.5 and 5.4, respectively) than the other eruptions (ratios are less than 2.4). Consequently, the 310 January retained eruptions are in the tropics and the July retained eruptions are shifted towards the 311 NH mid-latitudes (especially for 44 BCE) with the closer proximity of these eruptions to Greenland 312 balancing the reduction in poleward transport due to being in the summer hemisphere. The range in 313 retained SO2 emissions is also similar between seasons for both cases and therefore the similar 314 cumulative RF distributions despite differences in eruption latitude may be because of differences in 315 cumulative RF related to the position and strength of the tropical pipe (Fig. 5). 316 For all ten events, the injection height across all retained parameter combinations is not constrained, 317 although there are some combinations of SO2 emission and eruption latitude where injection height is 318 slightly constrained. Supplementary Fig. 4 shows that the January constrained parameter space is 319 often sloped, with lower emissions that have the highest injection heights, and higher emissions with 320 lower injection heights. The lack of constraint is important since although the time-integrated 321 deposition is generally not sensitive to the injection height, the cumulative RF is (Fig. 5). The 322 injection height is likely more important for the timing of the ice sheet deposition. 323 Our estimated mean values of plausible SO2 emissions for each eruption are generally higher than the 324 volcanic stratospheric sulfur injection (VSSI) estimates used to derive the PMIP4 prescribed sAOD 325 timeseries, EVA(2k) 18 (Table 1) values than the VSSI estimates, it is not surprising that the cumulative sAOD from EVA(2k) for each 331 of the ten eruptions is towards the lower end, or in the case of the 1230 CE eruption, almost outside of 332 our sAOD ranges. Our results consequently suggest that the EVA(2k) sAOD may be an 333 underestimate. Our uncertainty ranges are also generally higher (Fig 3b). 334 Because the natural meteorological variability (accounted for as a noise variance term on the 335 emulators) can cause further uncertainty on the predictions and can increase the probability that 336 parameter combinations are retained (for example at higher latitudes), constrained parameter 337 combinations obtained using only the emulator mean prediction of the deposition are also shown in 338 the supplementary information ( Supplementary Figs 6-9, Supplementary Table 1). Fewer plausible 339 eruptions are retained for each of the ten events, but the overall uncertainty on cumulative sAOD and 340 cumulative RF remains similar. 341 Table 1 Constrained ranges in SO2 emission (Tg SO2), eruption latitude (°N) and cumulative RF (MJ m -2 ) for the ten largest bipolar ice-core-derived sulfate 342 deposition signals in Greenland and Antarctica 4 (in rank order of magnitude). Year (BCE/CE) is the onset of the deposition signal. Also included are the 343 eVolv2k 18 Volcanic Stratospheric Sulfur Injections (VSSI) (plus or minus 2 standard deviations) (Tg of SO2). For all cases the minimum and maximum 344 retained injection plumes were 15-18 km and 25-28 km. Grey shading marks signals that are unconstrained because they are outside of our parameter space. 345 Although the latitude of 1815 Mt. Tambora is known, we keep the whole range in constrained parameters as an example for if the eruption had not been 346

350
We have constrained the eruption source parameters and the cumulative sAOD and cumulative RF of 351 eruptions corresponding to the ten largest bipolar ice-sheet sulfate deposition fluxes over the past 352 2500 years, without relying on transfer functions and scaling factors derived from single eruptions. 353 Our results suggest that there are many eruption realisations that could be consistent with ice-core-354 derived sulfate deposition fluxes, thus estimates of volcanic radiative forcing are more uncertain than 355 previously thought. 356 The cumulative RF has an uncertainty of at least ~300 MJ m -2 for the ten historical eruptions over the 357 last 2500 years that we have analysed, and can be as high as 424 MJ m -2 . To put this uncertainty range 358 in context, the cumulative RF predicted for the 1991 eruption of Mt. Pinatubo has been estimated to 359 be between -133 and -229 MJ m -2 (37) and is -203 MJ m -2 in the IPCC AR5 1 volcanic forcing series 360 (integrated over 1991 to 1996). Our uncertainty range on the cumulative RF of these ten past 361 eruptions therefore equates to approximately 1.5 to 3 times the total RF of 1991 Mt. Pinatubo. The 362 1991 eruption of Mt. Pinatubo led to up to 0.5°C of global surface cooling 45 , so we estimate that the 363 global mean temperature response to past eruptions could be uncertain to within at least 1.5°C, 364 assuming that temperatures respond linearly to forcing 1 . The eruptions have a range in cumulative RF 365 because there are many combinations of eruption latitude, SO2 emission, emission height and season 366 that produce ice-sheet sulfate deposition values that are consistent with measured anomalies. These 367 variations in eruption source parameters affect the amount of aerosol that is formed, the particle size 368 distribution, the horizontal and vertical distribution of the aerosol, its lifetime and hence its radiative 369 effect 37 . 370 Our constrained SO2 emissions are at the higher end of previous estimates (eVolv2k 18  attributing climate variability on millennial timescales because different volcanic forcings will lead to 380 different climate responses. Only one realisation of EVA(2k) with the medium VSSI predictions will 381 be used in PMIP4, thus none of this uncertainty will be accounted for. Our results also suggest there 382 are many other latitudes other than 0°N (to which unidentified eruptions leading to bipolar deposition 383 signals are assigned) that could lead to the deposition signals, and more often than not these are south 384 of the equator, especially if the eruption occurred in January. Hence, the assumption that eruptions 385 occurred at 0°N may not be realistic. 386 Our estimated radiative forcing uncertainty ranges are likely to be lower bounds for several reasons. 387 Our SO2 emissions were capped at 100 Tg; we considered eruptions only in January and July and 388 during the easterly QBO phase; we considered only one standard deviation uncertainty on the 389 emulator-predicted deposition and on the ice-core-derived ice sheet estimates during constraint; and 390 we did not include the emulator uncertainty on the sAOD and RF predictions, which was found during 391 validation to be very small (see Methods). 392 The difficulty in quantifying volcanic deposition anomalies in the model simulations and the 393 uncertainty in the emulator predictions are limitations of this study. It is possible that some of the 394 retained parameter combinations, especially at high latitudes, could be false positives that arise due to 395 variations in the non-volcanic sulfate deposition or due to the uncertainty in the emulator predictions. 396 However, this is also applicable to the real world. Our results have revealed that many eruptions could 397 be missed in the ice core records, or the volcanic sulfate deposition could be overestimated because of 398 internal variability and false attribution. (weeks). An emulator maps a model output (e.g. total sulfate deposited on Greenland) to the input 455 parameters (here the SO2 emission, eruption latitude and injection height) and can be used to predict 456 that model output for any combination of the input parameters that was not explicitly simulated. By 457 sampling from an emulator thousands of times, a multi-dimensional response surface of the model 458 output can be generated across the input parameter space, based on only a small set of model 459

simulations. 460
We generated four Gaussian process emulators 54-56 of the simulated deposition: total sulfate deposited 461 on Greenland following eruptions in January and July, and total sulfate deposited on Antarctica 462 following eruptions in January and July. Each emulator was constructed using R 57 and the 463 DiceKriging package 58 . Following a Bayesian statistical approach, each model response is assumed a 464 priori to follow a Gaussian process, which is then updated with the information in the model data easterly QBO phase), the standard deviation of the Greenland deposition was 6.2 kg SO4 km -2 (RSD = 505 20%) and the standard deviation of Antarctica deposition was 1.6 kg SO4 km -2 (RSD = 8%). Given 506 that the simulation of Mt. Tambora was initialized with a 60 Tg SO2 injection at the equator, it is 507 reasonable that an average noise variance term for our ensemble with emissions spanning 10 to 100 508 Tg SO2 and across both high and low latitudes, is higher. Regardless of the value of estimated noise, 509 we found that the overall shape and pattern of the emulated surfaces remained the same but the 510 emulator validation was improved using the given values. We found that higher estimates of this noise 511 variance led to poorer emulator validation and response surfaces with reduced variation in model 512 output versus the eruption source parameters. These values remain a best-estimate and reflect the 513 average variability in deposition across the whole parameter space and across both seasons. Although 514 it is possible to specify a heterogeneous noise term, it is likely this would introduce greater 515 uncertainty given a lack of information as to how internal variability terms may vary across the 516 parameter space. Because our simulations use prescribed SSTs, internal variability associated with 517 ENSO is not included. Similarly, we do not investigate variability associated with different QBO 518

phases. 519
We built four further Gaussian process emulators of the cumulative RF and cumulative sAOD for 520 each season (for eruptions occurring in January and July). These emulators were built without noise 521 because the sAOD and radiative forcing signals have relatively low variability compared to the 522 deposition (they are not determined by tropospheric meteorology) and validation of the emulators was 523 reasonable without an additional noise term. 524 Validation of the emulators is shown in Supplementary Figs. 11 and 12. The emulator predictions 525 follow the 1:1 line (marking a perfect prediction by the emulator) in all cases. The 95% confidence 526 bounds on the emulator predictions are larger for the deposition emulators ( Supplementary Fig. 11) 527 compared to the cumulative sAOD and RF emulators ( Supplementary Fig. 12) because the fit is more 528 uncertain and because of the additional noise term in the build. Overall, the emulators are reasonable 529 surrogates of the UM-UKCA output. Emulated response surfaces of the model outputs were produced 530 by sampling the predicted response of each emulator 1 million times over a three-dimensional grid 531 generated with 100 values of each eruption source parameter (covering the range in values simulated 532 for each parameter). 533

Constraining the eruption-realisations 534
Eruption-realisations are retained if the emulator-mean prediction of the Antarctica deposition plus or 535 minus one standard deviation and the emulator-mean prediction of the Greenland deposition plus or 536 minus one standard deviation overlaps with the ice-core-derived estimates and their uncertainty 4 . By 537 including the emulator uncertainty on the deposition emulators (which also accounts for ensemble 538 spread) and the uncertainty on the ice sheet composite observations, we provide a conservative 539 estimate of the range in eruption source parameter combinations and subsequently the range in 540 cumulative RF for the parameter space we have investigated. We repeat the constraint procedure 541 using the emulator-mean predictions only (without using the emulator standard deviation, so only 542 retaining combinations where the mean prediction lies within the observed range) to show the more 543 constrained estimates of plausible parameter combinations ( Supplementary Figs 6-9, Supplementary 544 Table 1). Most (except 44 BCE and 266 CE) of the ten ice-sheet deposition constraints we consider 545 are large enough in magnitude that it is unlikely that these signals could be produced by non-volcanic 546 sulfate anomalies in our simulations (i.e. they are much greater than 20 kg SO4 km -2 on Greenland and 547 10 kg SO4 km -2 on Antarctica where a clearer volcanic signal is observed in our simulations (Fig. 1)). 548 We constrain preindustrial sulfate deposition signals from simulations conducted using a present-day 549 atmosphere, where the background non-volcanic sulfate emissions are higher and large-scale 550 atmospheric circulation is faster. However, we find that the emulators predict the deposition following 551 the 1815 eruption of Mt. Tambora from preindustrial UM-UKCA simulations 40 and therefore this 552 difference is unlikely to significantly impact our results. 553 The EVA(2k) sAOD reconstruction is calculated using the Easy Volcanic Aerosol (EVA) forcing 554 generator 35 and Volcanic Stratospheric Sulfur Injections (VSSI) from the eVolv2k reconstruction 18 . 555 The sAOD timeseries for each of the 10 eruptions from the EVA(2k) reconstruction are shown in 556 Supplementary Fig. 13. We used five runs of EVA (with no background aerosol included) that were 557 run using the lower-end, middle and upper-end VSSI SO2 estimates and isolated our 10 eruptions 558 (upper and lower injections were calculated by summing/subtracting the one and two standard 559 deviation uncertainties from the middle best-estimate predictions). For each eruption we summed the 560 sAOD over 38 months to compare directly to the cumulative sAOD derived from the UM-UKCA 561 simulations. The resulting cumulative sAOD from the three runs are shown by the vertical lines in 562 Fig. 3. 563

Data Availability 564
Model output is included in Supplementary Table 2.  565