Volcanoes and ENSO: a re-appraisal with the Last Millennium Reanalysis

ABSTRACT


17
The El Niño-Southern Oscillation (ENSO), the quasi-periodic alternation of warm and cold phases of the tropical Pacific ocean- 18 atmosphere system, is the leading source of global interannual climate variability 1 . ENSO influences weather conditions 19 not only in the tropical Pacific 2 , but also globally through atmospheric teleconnections 3 . Skillful prediction of the ENSO 20 cycle, including its phase and amplitude, is therefore key to the successful forecasting of worldwide meteorological and 21 oceanographic conditions at sub-seasonal to seasonal scales.

22
External forcing has the potential to affect ENSO variability 4 . In particular, explosive volcanism may inject large amounts 23 of sulfate aerosols into the atmosphere, abruptly reducing incoming shortwave radiation and affecting the subsequent global 24 ocean-atmosphere climate variability for several years 5 . A causal relationship between large eruptions and ENSO would be a 25 significant source of climate predictability on interannual scales, and would be important for evaluating climate models' sensi-26 tivity to volcanic forcing, as well as assessing the risk of geoengineering solar radiation management schemes that emulate the 27 stratospheric sulfate aerosol loading characteristic of large explosive eruptions 6 . The link between volcanism and ENSO has 28 been vigorously debated since it was first proposed 7 . Since then, several tree-ring based observational studies have found sup-29 port for an El Niño-like response in the year following large eruptions 8-12 and at least five mechanisms have been proposed to 30 account for this relationship: (i) the ocean dynamical thermostat (ODT) 4,13,14 , which states that the upwelled water in the east-31 ern Pacific (EP) makes the region less sensitive to radiative forcing than the western Pacific, and leads to nonuniform Pacific 32 SST response to uniform incoming solar radiation reductions after eruptions; (ii) the land-ocean temperature gradient (LOTG) 33 mechanism 15-17 , which states that the low thermal inertia of the land introduces a LOTG after eruptions, which affects the 34 Pacific zonal wind anomalies and hence the ocean temperature; (iii) the subtropical wind stress curl mechanism 18, 19 , which 35 states that the initial enhanced cooling in EP after eruptions leads to a negative (anticyclonic) subtropical wind stress curl, 36 coral record from the heart of the tropical Pacific 25, 26 suggests no uniform ENSO response to all eruptions over the last mil-48 lennium, even for the largest eruptions 27 . This is in line with results from recent modeling studies using large ensembles that 49 allow quantification of the influence of stochastic as well as deterministic elements 24 . Indeed, ENSO is thought to be affected 50 by multiple uncertain or poorly constrained factors, including the phase of the quasi-biennial oscillation (QBO) 28 , the forcing 51 magnitude, location, and season of the eruption 19, 24 , as well as pre-conditioning of the ENSO state (neutral, Central Pacific 52 El Niño, Eastern Pacific El Niño, or La Niña) 16 . Yet, observational studies on this track are hindered by the limited num-53 ber of well known eruption events, the temporal resolution of volcanic forcing reconstructions, and the spatial and temporal 54 availability of proxy records. 55 Paleoclimate records offer a longer period of observation but conflicting accounts: reconstructions based mostly on tree-56 ring proxies [8][9][10][11][12] , which experience ENSO through teleconnections, have been used to argue of an El Niño-like response within 57 a year of the eruption; in contrast, reconstructions using corals from the core ENSO region 26, 27, 29 -which provide a firsthand, 58 albeit discontinuous, account of ENSO variations -do not support this conclusion. 59 In this study, we re-appraise the potential links between volcanism and ENSO, by integrating the latest paleoclimate 60 evidence from both tree rings and corals into a consistent dynamical framework, the Last Millennium Reanalysis (LMR) 30, 31 , 61 and interpret the results in the context of recent modeling work showing a large role for both initial and boundary conditions 62 in shaping the climate response to volcanism 16, 24 . 63 Corals vs tree rings 64 Coral archives are a natural choice for ENSO reconstruction due to their geographical proximity to ENSO centers of action 32 65 and the demonstrated link between the geochemistry of their skeletons and ENSO conditions 33,34 . However, their limited time 66 span produces discontinuous records, which can only be pieced together by splicing 26 . The longest and most complete coral 67 ENSO record published to date 26,27,35 is located at Palmyra atoll, at the edge of the Niño 3.4 region. The record covers 535 68 of the past 1000 years, leaving many gaps, and with limited replication over common intervals. Tree-ring based proxies, on 69 the other hand, have been used to build long, cross-dated, and heavily replicated reconstructions that continuously span the 70 Common Era. Their distance to ENSO centers of action means that they rely on teleconnections between the tropical Pacific 71 and their local terrestrial rainfall and temperature anomalies, leaving them vulnerable to confounding factors. Combining 72 data from both archives could ameliorate their individual limitations. Before doing so, however, we compare reconstructions 73 assimilating the two archives separately, as this provide insights into the possible causes of the discrepancy between previous 74 studies and on how to interpret a combined signal.

75
Our strategy leverages the data assimilation algorithm of the Last Millennium Reanalysis (Methods). We first assimilate 76 corals to reconstruct tropical Pacific surface temperature during the boreal winter (DJF).

77
Our coral collection includes the synthesis of Ocean2k 29, 36 , supplemented by the latest Palmyra record 27 . The resultant 78 reconstruction is denoted as LMR (Corals), and its spatial and temporal skill is presented in Fig. 1(a,d). The skill of the 79 Niño 3.4 reconstruction is remarkably high compared to many other existing reconstructions (see Extended Data Fig. 1), with 80 a temporal correlation coefficient (R) of 0.86 (that is, 74% of shared variance) and a coefficient of efficiency (CE) 37 of 0.71 81 against ERSSTv5 (Extended Reconstructed Sea Surface Temperature v5) 38 . However, due to the temporal gaps in the Palmyra 82 record, the reconstruction is incomplete over the last millennium (Fig. 2b) and the ensemble-mean variability collapses when 83 Palmyra observations are unavailable.

84
The assimilation of tree rings for Niño 3.4 reconstruction is challenging due to their distance to the target region. To 85 overcome this, we gather the six best tree-ring based Niño 3.4 predictors identified by Li et al. 10  in year 1 after instrumental era eruptions (as the colored dots tend to hover above the linear model fit), yet this behavior is not 92 obvious in the reconstruction LMR (Li13b6) (Fig. 1h), suggesting an advantage for processing the proxy observations within 93 the LMR paleo-data assimilation framework. The spatial verification (Fig. 1b) indicates that the reconstruction skill peaks 94 around the center of the Niño 3.4 region and decays quickly away from it. A major advantage of the tree-based reconstruction 95 is its continuous nature over the past millennium (Fig. 2c). 96 We next assess the ENSO response to volcanism in both reconstructions. To perform an equitable comparison, we select 97 large eruptions defined as volcanic stratospheric sulfur injection (VSSI) greater than 6 Tg S according to the eVolv2k version 98 3 dataset 43 , and focus on periods when the Palmyra coral record is available (Fig. 2a). We assess the forced ENSO responses 99 using both the widely-used superposed epoch analysis (SEA) approach (also known as compositing), as well as a ranking 100 analysis (see Methods). The SEA (Fig. 3a,b) applied to the 12 large events when Palmyra is available (Fig. 2a) suggests that 101 there is no Niño 3.4 composite value that is significantly higher than that of randomly drawn non-volcanic years, even at the 102 relatively permissive 90% level.

103
However, SEA is sensitive to dating uncertainties, as the compositing procedure requires that the reconstruction segments 104 of individual events align precisely, without which otherwise minor timing offsets can compound and damp the composite 105 signal. SEA may also be affected by the small sample size. This may first come into play through a type I error, which is 106 tested in the SI Text 2. We find that even small ensembles (n = 5) can mitigate this problem, so it is not material for this 107 particular issue. However, like all averaging methods, an SEA carried over a small sample may be dominated by a small 108 number of events with extreme anomalies, with most events showing a modest response 44 . Thus, a "significant" composite 109 response is a necessary but insufficient condition to establish a physical response. Only when most of the events show strong 110 responses can we be confident about the robustness of the relationship. Therefore, it is important to parse composites into 111 the contributions from, and consistency of, individual events. This is achieved through a ranking analysis ( Fig. 3c-f), which 112 compares the Niño 3.4 response of each event to the distribution of all non-volcanic years. Events with response values larger 113 than a certain quantile (say, 80, 90, or 95%) of the distribution of non-volcanic years are considered "significant".

114
Though the ranks of eruption-year events in the two reconstructions do not agree with each other in detail, both LMR 115 (Corals) and LMR (Li13b6) point to the same conclusion: that most eruption years are not statistically different from non-116 eruption years in both year 0 and year 1, at any significance level.

117
Combining corals and tree rings 118 Given this agreement, we combine both archives to perform a continuous reconstruction, denoted as LMR (Corals+Li13b6).

119
Verification statistics (Fig. 1c,f) show a slight skill improvement compared to LMR (Corals), with R = 0.86 and CE = 0.71 120 against ERSSTv5. More importantly, the temporal gaps in LMR (Corals) over the past millennium are now filled with the 121 information from Li13b6 (Fig. 2), which allows a larger sample size of eruption events, hence a more robust evaluation of 122 ENSO's response to volcanism.

3/20
To assess individual events, we compare rankings using all eruptions vs pre-instrumental eruptions only (Extended Data 143 Fig. 3c,d). It can be seen that the instrumental period eruptions overall show a stronger Niño 3.4 anomaly than the pre-144 instrumental period eruptions (Extended Data Fig. 3c). Without instrumental-era eruptions, only (3, 2, 1) out of 15 events 145 show a larger Niño 3.4 response than (80%, 90%, 95%) of the non-volcanic years (Extended Data Fig. 3d), indicating that 146 most of the pre-instrumental period eruptions are not significantly different from non-volcanic years.

147
The above analysis suggests that the reconstructed ENSO response to volcanism is non-stationary, at least between pre-148 instrumental and instrumental period eruptions.

149
This could be caused by data attrition, which reduces variability of a reconstruction going back in time since the proxy 150 observations provide the only source of variability in our reconstructions. Moreover, MADA is compiled using a correlation-151 weighted, ensemble-based modification of the "point-by-point regression" 40 , and the ensemble members become less similar 152 to each other during earlier periods, over which temporal variability is damped compared to recent intervals. In addition, 153 dating accuracy is a non-negligible problem for fossil corals 45 , where errors compound back in time, potentially reducing 154 the variability of composite series derived from them 34 . In addition, we note that the instrumental period is quite short and 155 devoid of the very large eruptions that occur in the pre-instrumental period, so results obtained over this period may not be 156 statistically representative. Furthermore, instrumental records have shown evidence of the potential of coincidence between 157 ENSO activity and volcanism, at least for the Agung (1963), El Chichon (1982, and Pinatubo (1991) eruptions, when strong 158 El Niño events were already underway before the eruptions went off 46 .  Effects of forcing asymmetry 164 Previous studies have suggested that the hemispheric asymmetry of volcanic forcing is another factor that may differentially 165 affect ENSO activity 15,20,21,24 . Therefore, we investigate the relationship between forcing parameters and the Niño 3.4 SST 166 from LMR (Corals+Li13b6). Tropical eruption events are categorized by different hemispheric asymmetry levels (below 167 0.8, between 0.8 and 1.5, and above 1.5) defined as aerosol spread based on ratio of Greenland to Antarctic sulfate flux 43 .

168
Extreme events are defined as having a VSSI>20 (1230,1257,1458,1815). Extended Data Fig. 4   Using state-of-the-art datasets and methods, combining the strengths of both coral and tree-ring records, the present evidence 177 is as yet unable to detect an effect of explosive volcanism on ENSO phase, consistent with D20 27 . This conclusion is supported 178 by corals and trees independently. Yet, corals and trees disagree regarding the relative rank between individual events ( Fig.   179 3). This issue might be caused by attrition in both networks, which reduces reconstruction variance back in time. Similarly, 180 temporal gaps in the Palmyra coral record limit the pool of eruption events and non-volcanic years to a smaller size, making 181 the comparison less stable.

182
Beyond the issue of attrition, proxy archives harbor limitations of their own. The tree-ring records used here, for instance, 183 are sensitive to the local hydrological expression of volcanic eruptions, which may be mistaken for the remote effects of ENSO 184 via teleconnections 22 . The ENSO reconstruction based on the South American Altiplano composite 42 also suggests that the 185 non-stationary behavior of ENSO variance is masked by intrinsic tree-ring width variability. Indeed, when ENSO variance is 186 low, the teleconnection is weak, and tree ring variability is more reflective of local temperature and moisture conditions than 187 those from the remote tropical Pacific; such non-stationarity may obscure the volcanoes-ENSO link. The current coral network 188 is more proximal to ENSO centers of action, yet is dominated by records from the western and central equatorial Pacific (Fig. 189 1a) which capture La Niña events more faithfully than El Niño events, and tend to underestimate the amplitude of large eastern RSST highlights the impact of volcanism on ENSO relative to the tropical mean cooling, and is computable in our recon-195 structed SST fields. This is shown in Extended Data Fig. 5, which displays a temporally-flat tropical average temperature 196 anomaly for most of the years between 1100-2000 CE, resulting in RSST-based Niño 3.4 anomalies that are indistinguishable 197 from the SST-based ones. This is a direct result of the sparse coral network before the nineteenth century, and will not im-198 prove until this network is vastly expanded over tropical oceans. However, even with the RSST conversion, we still observe 199 consistently low significance ratio in the ranking analysis of the last millennium climate model simulations (see SI Text S4 for 200 details).

201
Another caveat comes from limitations associated with the choice of the model prior in the LMR paleo-DA framework.

202
The model prior here refers to the model simulation from which climate states are chosen at random by the Kalman filter 203 algorithm. Its chief role is to provide the spatial covariance information within and between fields, and affects how the 204 information from proxy observations propagates to locations where those observations are not available. According to a 205 recent study 49 , the known biases in the location of the South Pacific Convergence Zone (SPCZ) in most climate models leads 206 to incorrect inferences about Niño 3.4 SST from corals located in the SPCZ region, in such a paleoclimate data assimilation 207 context. Moreover, ref. 49 shows that corals located in both the SPCZ and Niño 3.4 regions produce local cooling during 208 the 1809 and 1815 eruptions, but all prior ensembles considered (including those drawn from 20th century reanalyses) have a 209 covariance pattern that yields remote influence (that is, the influence of one region on another) inconsistent with this pattern.

210
In addition to model bias, this suggests that model priors conditional on eruption time would be important for properly 211 representing the information in proxy records. We note that, in the present study, these biases are mitigated prior to 1800, 212 when only proxies tied to the Niño 3.4 index are used.

213
Finally, a recent modeling study 24 suggests that the ENSO response to volcanism is rather weak during DJF because it 214 changes sign around that season. Yet, the response could be strong before the sign changes, which is usually during January-215 September (JAS) and/or October-December (OND). Repeating our analysis using JAS and OND as target seasons, we find 216 our conclusions insensitive to this choice (SI Text 5). This is in apparent contradiction to that modeling study 24 , yet may 217 be explained by the limitations of currently available proxy records noted above. Assigning eruption events to specific years 218 based on volcanic forcing reconstructions could also lead to time offsets in analyses using different target seasons 50 , which 219 will potentially obscure the relation to eruption events. This is a problem that we cannot resolve currently without a better 220 knowledge of the eruption season.

221
That we do not find a consistent ENSO response in the relatively small statistical sample of the last millennium is not 222 evidence that there is zero predictability associated with volcanoes. As suggested by recent modeling studies 16, 24 , the forcing 223 magnitude, location, and season of the eruption, as well as pre-conditioning of the ENSO state can greatly affect the ENSO 224 response to volcanic eruptions. These multiple factors must align to favor the development of ENSO events, providing a source 225 of predictability only when these factors are known with a sufficient degree of accuracy. It is therefore reasonable to expect 226 that the lack of an observed, consistent relationship between ENSO phases and explosive volcanism in our reconstructions 227 may be due, in part, to an imperfect knowledge of these factors. However, even when controlling for eruption timing in PMIP3 228 and CESM-LME simulations, a ranking analysis still demonstrates an inconsistent ENSO response to volcanism when this 229 major source of variability is held fixed (SI Text 4), suggesting that all factors need to be jointly determined, or that a larger 230 ensemble is needed to discern common trends.

231
Given the large number of degrees of freedom, a large sample size is needed to isolate a consistent signal -larger perhaps 232 than offered by the past millennium. It is thus important to develop more high-resolution proxy records spanning the tropical 233 oceans over the last millennium, and possibly extend them through the longer Holocene. The reconstruction of volcanic 234 forcing also needs to be expanded with longer temporal coverage for more robust statistics. Until these goals have been 235 achieved, it is unclear how much one can conclude from the paleoclimate record about the contribution of volcanic eruptions 236 to ENSO dynamics, or to the assessment of the risks posed by solar radiation management strategies in relation to ENSO.             covariance localization 56 radius of 25,000 km is applied, and we note that the results are insensitive to this exact choice.

365
Note that the calibration period for the PSMs is 1850-2000 CE so as to achieve the best reconstruction skill. To guard 366 against potential overfitting, as well as the potential impact of climate change we see in NADA PC2 and the Kauri composite 367 (Extended Data Fig. 2d,e), a cross-validation of the reconstructed Niño 3.4 is performed with disjoint calibration and validation 368 periods (see SI Text 1); the result suggests that the reconstruction skill is stable to this choice.

369
Data Sources 370 We consider information from seasonally-sensitive or monthly-resolved proxy records. The other predictor is the west-central Argentina composite 57 , which contributes negligibly to reconstruction skill (with an 380 explained variance of Li13 of 1.8%), and is ignored here. We refer to these best six records as "Li13b6".

381
We reproduced the principal components of NADA and MADA via principal component analysis (  For each reconstruction, dark shading denotes the interquartile range, and light shading denotes the central 95% region, from 2.5% to 97.5%. R=correlation coefficient, CE=coefficient of efficiency 37 . (g-i) Scatter plot of the data points in (d-f). The grey dashed curve represents the linear regression fitting curve. The black and colored dots denote the data points at year 0 and year 1 relative to large eruption years (1883,1902,1913,1951,1963,1982,1991)  Signif. ratio: (4,2,1)/22 Year 1 f 80% 90% 95% 50% Fig. 4. Same as Fig. 3, but for LMR (Corals+Li13b6) regarding the 12 events when Palmyra is available and all the 22 events over the past millennium.    (n=22) excluding 1883, 1902, 1913, 1951, 1963, 1982, 1991 (n=15)   For each reconstruction, dark shading denotes the interquartile range, and light shading denotes the central 95% region, from 2.5% to 97.5%. R=correlation coefficient, CE=coefficient of efficiency37. (g-i) Scatter plot of the data points in (d-f). The grey dashed curve represents the linear regression fitting curve. The black and colored dots denote the data points at year 0 and year 1 relative to large eruption years (1883,1902,1913,1951,1963,1982,1991) as in Li et al.10, respectively. Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors. Figure 2 (a) The 22 large eruption events defined as the volcanic stratospheric sulfur injection (VSSI) greater than 6 according to eVolv2k version 343. 12 events over the years when the Palmyra coral record26, 27, 35 is available are colored in black, while other events are colored in grey. See Extended Data Table 1 for the details of the metadata. (b-d) Same as Fig. 1d-f, but for the past millennium. Figure 3 (a-b) Superposed epoch analysis (SEA) of the reconstructions LMR (Corals) and LMR (Li13b6) regarding the 12 events when Palmyra is available. Solid curves with dots denote the composite mean, and the light dots denote the Niño 3.4 anomaly at each year for each individual event. The light grey dashed curves denote the 1%, 5%, 10%, 90%, 95%, and 99% quantiles of the composite means from 1000 bootstrap draws from non-volcanic years (Methods). (c-d) Ranking analysis of the LMR reconstructed Year 0 Niño 3.4 values. The grey shaded area denotes the distribution of the Niño 3.4 anomaly value over all nonvolcanic years, whose 50%, 80%, 90%, and 95% quantiles are denoted by vertical dot-dashed curves, serving as significance levels (Methods). The vertical solid lines mark individual volcanic events; for each, the horizontal axis position denotes the Niño 3.4 anomaly value, and the vertical axis position denotes the relative rank of the Niño 3.4 anomaly value compared to all other events. The circle/downward triangle/upward triangle/diamond marker represents that a volcanic event has a Niño 3.4 anomaly value that is below 80%/between 80-90%/between 90-95%/above 95% of that over the non-volcanic years. The significance ratio denotes the number of events that are above the 80%, 90%, and 95% significance levels, respectively, out of all volcanic events. (e-f) Same as (c-d), but for the year 1 Niño 3.4 values.

Figure 4
Same as Fig. 3, but for LMR (Corals+Li13b6) regarding the 12 events when Palmyra is available and all the 22 events over the past millennium.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SI.pdf