A “La Niña-like” state occurring in the second year after large tropical volcanic eruptions during the past 1500 years

Using an ensemble of nine El Niño/Southern Oscillation (ENSO) reconstructed proxies and volcano eruption proxies for the past 1500 years, this study shows that a significant La Niña state emerges in the second year (year (2) hereafter) after large tropical volcanic eruptions. The reasons for the development of La Niña are investigated using the Community Earth System Model (CESM). In the volcanic eruption experiment (Vol), a robust La Niña signal occurs in year (2), resembling the proxy records. The eastward positioning of the western North Pacific anomalous anticyclone (WNPAC) in Vol plays a critical role in the advanced decay of year (2) warming and the strong intensification of cooling in the equatorial eastern Pacific. The enhanced easterlies located on the southern edge of the WNPAC can stimulate consecutive oceanic upwelling Kelvin waves, shallowing the thermocline in the eastern Pacific, thereby resulting in a greater cooling rate by the enhanced thermocline feedback and cold zonal advection. Over the equatorial eastern Pacific, the reduced shortwave radiation contributes to the advanced decay of warming, while the upward latent heat flux augments the strong intensification of the cooling. Essentially, the eastward positioning of the WNPAC is a result of the volcanic forcing. The volcanic effect cools the maritime continent more than its adjacent oceans, thus pushing convective anomalies eastward during year (1). This induces vertical thermal advection and upward surface latent heat flux, thereby suppressing the development of warm Sea Surface Temperature over the central-western Pacific and causing the eastward positioning of the WNPAC in Vol.


Introduction
The El Niño/Southern Oscillation (ENSO) plays a dominant role in global interannual climate variability (Webster et al. 1998) through atmospheric teleconnections, regulating extreme climate events which have profound social and economic impacts (Wang et al. 2000(Wang et al. , 2013. While ENSO is primarily determined by atmosphere-ocean interactions within the coupled climate system (Bjerknes 1969;Philander et al. 1984;Cane et al. 1986), recent studies of long-term paleoclimate records suggest that large tropical volcanic eruptions can induce an El Niño-like state, or increase the occurrence of El Niño events (Adams et al. 2003;Mann et al. 2005;McGregor et al. 2010).
The process by which tropical volcanic eruptions influence ENSO has been studied via numerical simulation (Mann et al. 2005;Emile-Geay et al. 2008;Zanchettin et al. 2012;Ohba et al. 2013;Ning et al. 2017;Predybaylo et al. 2017;Stevenson et al. 2017). Mann et al. (2005) explained Sea Surface temperature (SST) change over the eastern Pacific after volcanic eruptions by utilizing an air-sea coupled model (Zebiak and Cane 1987), invoking the dynamical thermostat hypothesis (Clement et al. 1996). This hypothesis argued that cooling occurrs more rapidly in the western Pacific than the eastern Pacific with the decreasing of surface shortwave radiation in the tropical Pacific, resulting in reduced zonal SST gradients. Meanwhile, Ohba et al. (2013) explained the occurrence of El Niño by a land-sea cooling contrast hypothesis. They demonstrated that the rapid reduction of surface temperature over the maritime continent caused by the large volcanic eruptions may weaken the Walker circulation, leading to decreased subsurface water upwelling and then contributing to a warm phase of ENSO.
Taking the aforementioned into consideration, a central question for exploration is: What happens to the Pacific SST after initial El Niño events are induced by large tropical volcanic eruptions? Adams et al. (2003) proposed that the El Niño occurs after the tropical eruptions, followed by a La Niña-like state. Their reconstructed results suggested a significant 'rebound' from an El Niño into a La Niña-like state. However, the subsequent La Niña-like rebound seems to be weaker than the initial El Niño response. Recently, significant La Niña-like rebounds have been found in the instrumental SST composites around nine strong volcanic eruptions since 1870 (McGregor and Timmermann 2011). This result seems to be consistent with other studies (e.g., McGregor et al. 2010;Wilson et al. 2010). Further, Maher et al. (2015) found a La Niña-like response two years after the tropical volcanic eruptions in 122 historical ensemble members from the Coupled Model Intercomparison Project 5 (CMIP5).
It should be noted that the record durations documenting the second-year La Niña are relatively short, and do not contain many cases of tropical volcanic eruptions. Furthermore, some proxy evidence indicates that the anomalous eastern Pacific warming in the first year of volcanic eruptions is followed by a neutral SST state after year 2 (Li et al. 2013). Moreover, uncertainties exist in different proxy data and the reasons for divergent conclusions need to be further explored. Therefore, it is necessary to gather more ENSO proxy data on a longer time-scale to establish if this La Niñalike rebound exists. Meanwhile, the specific processes and mechanisms in action during the La Niña development phase under volcanic forcing remain to be revealed.
The positive feedback explaining the growth of El Niño events cannot explain the ENSO phase transition from a warm to a cold phase (Bjerknes 1969). The turnabout of El Niño to La Niña has been explained by delayed oscillator theory (e.g. Suarez and Schopf 1988;Battisti and Hirst 1989;Cane et al. 1986). The key process of this theory is the delayed effect induced by the upwelling Rossby waves excited in the central Pacific during El Niño, which propagate westward and generate upwelling Kelvin waves upon reflection from the western boundary of the Pacific. These upwelling Kelvin waves further adjust the thermocline and finally shut down the warming in the eastern Pacific. The other hypothesis involves a discharge-recharge paradigm in the zonal mean thermocline variation (Jin 1997;Li 1997). Recently, instrumental records have suggested that the rapid warm-to-cold transition phase of ENSO is attributed to both dynamic (atmospheric wind responses in the western Pacific) and thermodynamic (heat flux) forcing (Chen et al. 2016). The former is related to the western North Pacific anomalous anticyclone (WNPAC) during the mature phase of El Niño, which is maintained by an interaction between off-equatorial atmospheric Rossby waves and the underlying anomalous SST (Wang et al. 2000). The Philippine Sea anticyclonic anomaly deepens the thermocline locally, and the associated equatorial easterly anomalies trigger upwelling Kelvin waves, raising the thermocline and enhancing cooling in the equatorial eastern Pacific (Wang et al. 1999). The latter involves the cloud-radiation-SST feedback and evaporation-SST feedback caused by SST forcing-induced the nonlinear atmospheric wind, convection, and moisture (Chen et al. 2016). This rapid ENSO transition is mainly controlled by internal variability in the climate system. Liu et al. (2017) suggested a rapid warm-to-cold transition would occur in the second year (year (2) hereafter) after volcanic eruptions due to the easterly anomalies over the western Pacific using Community Earth System Model (CESM), but the underlying mechanism is still unclear.
In this study, nine sets of conventional cold-season ENSO proxy data are collected and examined, demonstrating the existence of a transition from El Niño to La Niñalike rebound after tropical volcanic eruptions. Particularly, it seeks to understand the specific processes during this transition from El Niño to La Niña by analyzing the past 1500 years of numerical simulation results with CESM. This paper is organized as follows. Section 2 introduces the model, experimental design, and proxy data. Section 3 presents the temporal and spatial evolution of SST over the eastern Pacific after large tropical volcanic eruptions. Section 4 offers an investigation of the physical processes controlling La Niña's development in the year after the volcanic eruptions. Section 5 presents the volcanic influence on the eastward shift of the location of the WNPAC, a key system causing the La Niña rebound. Finally, Sect. 6 summarizes the conclusions.

Model and experiments
The model used in this study is the Community Earth System Model (CESM) with the Parallel Ocean Program version 2 (POP2) and the Community Atmosphere Model, version 4 (CAM4) (Neale et al. 2013), developed by the National Center for Atmospheric Research in July 2010. The low-resolution version (T31_g37, which is equivalent to 3.75° × 3.75° in the atmosphere and land components) has been used to carry out simulations for the past 1500 years because of the long time required to integrate the model and the limited computing resources at our disposal. This means that CAM4 has a global range of 48 grid points in latitude, 96 grid points in longitude, and 26 vertical levels, while POP2 has 116 grid points in latitude, 100 grid points in longitude, and 60 vertical levels. The CESM is the one of the main models cited in the Fifth Assessment Report of the Intergovernmental Panel on Climate Change and the performance of CESM has been verified by many researchers Ning and Bradley 2016;Otto-Bliesner et al. 2016;Yan et al. 2016). CESM is used to simulate a long-time integration, including a control experiment (Ctrl) and a volcanic sensitivity experiment (Vol). The Ctrl, which consists of a pre-industrial 400-year spin-up run and a 2000-year control simulation, is run by the fixed external forcing condition of 1850 A.D (Wang et al. 2015). The Vol run, which spans from 500 A.D. to 2000 A.D., is based on the year 2400 in Ctrl under the initial forcing of 1850. The volcanic forcing used to drive the model is based on Ice-core Volcanic Index 2 (Gao et al. 2008), in which the volcanic forcing is a function of month, latitude and altitude. The volcanic forcing is the only changing external forcing during the 1500-year simulation in Vol.

Proxy data
The proxy data for ENSO used in this research are listed in Table 1. The first three reconstructed data in Table 1 are derived from tree ring records in North America, while the other six are the proxies from mixed sources, such as tree rings, corals and ice cores. There are five records before 1600 and the longest one is from 900 to 2002 (Li et al. 2011). All of these ENSO records are highly correlated with the instrumental SST over Niño 3.0 region (5°S-5°N, 150°-90°W) in HadISST during the industrial period (Rayner et al. 2003) (at the 99% confidence level).
These nine ENSO reconstructions are used to detect the response of ENSO to the strong tropical volcanic eruptions. To date the year of tropical volcanic eruptions, the latest records of atmospheric aerosol loading (Sigl et al. 2015) are used and the events of large tropical volcanic eruptions (with annual-mean global aerosol amount above 4 Tg in the eruption year) in Vol are extracted, following the classification defined by Liu et al. (2016).

Method
The superposed epoch analysis (SEA) (Haurwitz and Brier 1981) was used to identify the volcanic influence on warmto-cold ENSO transition and the related atmosphere-ocean interaction. Each window includes 7 years before and 8 years after each eruption. Following Adams et al. (2003) and Liu  (2016), in order to avoid the possibility that any single volcanic eruption may dominate the epochal signal, the data were normalized by the magnitude of the largest absolute anomaly magnitude in each window. Monte Carlo resampling (ensembles of 10,000 random numbers from the composite matrices) was used to detect significant changes (Adams et al. 2003).
In this study, a budget analysis of the tendency of the mixed layer temperature anomalies (MLTA) is performed to examine the physical processes in the formation of La Niña. The MLTA tendency equation can be derived as follows: where the overbars and primes indicate climatological and anomalous quantities respectively; T denotes the mixed layer temperature. As the tropical volcanic eruptions can induce tropical surface cooling, the mean tropical (30°S-30°N) mixed layer temperature anomalies in T are removed. Further, u, v, and w represent the zonal and meridional currents, and upwelling velocities respectively; ∕ x , ∕ y , and ∕ z represent the 3D gradient operator; Q net denotes the net downward heat flux at the ocean surface; (= 10 3 kg m −3 ) is water density; C p (= 4000 J kg K −1 ) is the specific heat of water; H represents the mixed layer depth, here taken as the top 70 m; and R denotes the residual term.

Observed and simulated SST changes in the equatorial Pacific
The superposed epoch analysis (SEA) results are shown in Fig. 1. The confidence intervals of SEA are calculated by a total of 10,000 Monte Carlo simulations. Year '0' represents the immediate cold season after volcanic eruptions and year '− 1' shows the initial ENSO state before tropical eruptions. Seven out of 9 proxy reconstructions reveal a significantly decreased SST over the eastern Pacific in the cold season of year (2) after the large tropical volcanic eruptions, representing a strong La Niña or La Niña-like state. Proxy data reconstructed by Li et al. (2013) (Fig. 1e) show a significant El Niño-like response in year (1), followed by a neutral state in years (2) and (3), not a La Niña-like state. The uncertainties in different reconstruction methods may be one of the reasons for this mismatch. In order to reduce these uncertainties, the ensemble of the nine cold-season ENSO reconstructions is used. The ensemble results (Fig. 1j) clearly show a transition from the year (1) El Niño-like state to the year (2) La Niña-like state after the eruptions. Notwithstanding, although the proxy data reveal the La Niña-like response in year (2), it remains unclear what caused the La Niña rebound. Thus, the model results are used to elucidate the specific processes and mechanisms responsible for the warm-to-cold ENSO transition after tropical eruptions. Firstly, the model's performance needs verification. The transition from the year (1) El Niño state to the year (2) La Niña state is also well captured by the SEA results over the Niño 3.0 region (Fig. 1k) and the Niño 3.4 region (5°S-5°N, 170°-120°W) (Fig. 2a) in the Vol experiment. Although the decreased SST which occurred over the eastern Pacific the immediate winter in the year of volcanic eruption does not match the reconstructed counterpart. This disagreement also appears in previous model simulations (e.g., McGregor and Timmermann 2011;Stevenson et al. 2016). The simulation shows that the tropical Pacific SST decreases rapidly in the immediate cold season after the large tropical eruptions, which may lead to cooling over the eastern Pacific ( Fig. 1k) but there remains an El Niño-like SST gradient across the Pacific.
To extract the volcanic-induced ENSO transition changes, both the Vol and Ctrl experiments are needed for the following analysis. The temporal evolutions of the SST anomalies over the Niño 3.4 region (5°S-5°N, 170°-120°W) are shown in Fig. 2. The composite El Niño episodes decay rapidly and turn into a La Niña phase in year (2) after the volcanic eruptions (Fig. 2a). This composite evolution shows a strong La Niña with a minimum of about − 0.7 °C. However, it is noted that the volcanic aerosol can cause an overall cooling in the tropics. In order to detect the volcanic influences on the ENSO transition, the tropical mean cooling caused by the volcanic eruption was removed. Here the negative SST anomalies averaged over the tropics (30°S-30°N) were used (indicated by the blue curve in Fig. 2a) as the volcanicinduced mean cooling state. As the composite mean SST internal variation over the Niño-3.4 region is about 0.87 °C at El Niño's peak in the Vol experiment (hereafter Vol for short), the El Niño year in the Ctrl experiment (hereafter Ctrl for short) that corresponds to its amplitude in Vol was selected, resulting 73 El Niño events in Ctrl. Figure 2b compares evolutions of the Niño-3.4 SST anomalies simulated in Ctrl and Vol. El Niño's decay occurs about two months earlier in Vol than in Ctrl (starting from Dec (1)). The decreasing tendency of the anomalous SST during the decay phase of El Niño (Dec (1) to Mar (2)) is -0.19 °C month −1 in Vol, while − 0.05 °C month −1 in Ctrl. There is a significant difference between them according to the Mann-Kendall test (at the 95% confidence level) (Kendall 1955;Wang and Ding 2006). Additionally, in the later part of year (2), the intensity of La Niña in Vol becomes stronger than that in Ctrl. The decreasing tendency of the anomalous SST during the developing phase of La Niña (Sep (2) to Dec (2)) is − 0.08 °C month −1 in Vol, while − 0.03 °C month −1 in Ctrl. There is also a significant difference between them after the Mann-Kendall test (at the 95% confidence level). This means that the significant advanced decay of warming occurs at the beginning of year (2), and the strong re-intensification of the cooling happens in the later part of year (2) over the equatorial eastern Pacific. Figure 3 compares the evolutions of the horizontal distribution of SST anomalies associated with the intrinsic ENSO mode in Ctrl and the ENSO mode induced by volcanic eruptions in Vol. Similar to the intrinsic El Niño mature phase in Ctrl (Fig. 3a), the warming in Vol matures over the eastern equatorial Pacific in the winter of year (1) (Fig. 3b). This warming event after volcanic eruptions has been reflected in many proxy data (Fig. 1) and previous simulations (Ohba et al. 2013;Maher et al. 2015;Stevenson et al. 2016). However, the location of the maximum warming in the El Niño mature phase in Vol (Fig. 3b) is located about 40° longitude to the east of that in Ctrl (Fig. 3a). At the same time, the cooling center over the WNP (0°-20°N, 120°E-160°E) in Vol is stronger and has a larger spatial extent than those in Ctrl.
In year (2), El Niño decays rapidly under the volcanic forcing. By spring (2) warming over the central equatorial Pacific is very weak (Fig. 3d). In summer (2), the negative SST anomalies occur across the entire equatorial Pacific. Thereafter, the cooling phase develops over the equatorial central-eastern Pacific (CEP) and reaches its mature La Niña state in winter (2/3) (Fig. 3j). The warming over the eastern Pacific decays earlier in spring (2), and the amplitude of the cooling in winter (2/3) is larger than that in Ctrl (Fig. 3i).
Essentially, such findings lead to the central questions: Why does El Niño decay earlier in Vol than in Ctrl? And,  (Stahle et al. 1998;Mann et al. 2000;D'Arrigo et al. 2005;Braganza et al. 2009;Cook et al. 2008;Li et al. 2011Li et al. , 2013McGregor et al. 2010;Wilson et al. 2010), j is the ensemble of all these nine reconstruc-tions, and k results from the CESM simulation. Year '0' represents the immediate cold season following the volcanic eruptions. The red color bars mark the first and second year after the eruptions. Confidence limits (90, 95%) based on Monte Carlo resampling (n = 10,000) are marked by the horizontal lines why does La Niña reach a stronger intensity in Vol than in Ctrl?

Physical processes responsible for the La Niña development
This section attempts to address the questions previously raised.

Oceanic processes determining La Niña's development
The discharge-recharge process can be an important mechanism leading to the ENSO transition (Jin 1997;Li 1997).
To examine the recharge-discharge paradigm during the warm-to-cold transition year, the heat content defined by the warm water volume (WWV) above the 18 °C isotherm between 5°N-5°S, 120°E-80°W was computed (Meinen and McPhaden 2000). Water above the 18 °C isotherm is defined as "warm water" in this study. As shown in Fig. 4, the WWV reaches the maximum before June (1) in both Vol and Ctrl. The discharge process in Vol leads that in Ctrl. This means that the earlier meridional divergent Sverdrup transport in Vol, which shallows the zonal mean thermocline depth, causes the advanced development of La Niña, consistent with the results in Figs. 2 and 3. Meanwhile, the discharge process is stronger in Vol than that in Ctrl from Sep (2) to Dec (2), conceivably causing the stronger La Niña.
To better understand the physical processes during the warm-to-cold transition after the volcanic eruptions, a budget analysis of MLTA averaged over the Niño 3.4 region, was completed. The relative contributions of each of the composite mixed layer terms during El Niño's decaying phase (Dec (1) to Mar (2)) are shown in Fig. 5. The MLTA tendency from Dec (1) to Mar (2) is − 0.05 °C month −1 in Ctrl and − 0.18 °C month −1 in Vol (Fig. 5a). The difference between them is significant after conducting a Mann-Kendall test (at the 95% confidence level). The vertical advection terms play a crucial role during the advanced decay of El Niño (Fig. 5b). In particular, the anomalous temperature gradient by mean upwelling ( −w T � ∕ z ), closely related to the thermocline feedback (An and Jin 2001), is − 0.10 °C month −1 under the volcanic forcing, stronger than that (− 0.02 °C month −1 ) in Ctrl. It means that the mixed layer can be shallower over the equatorial eastern Pacific in Vol. The zonal advection also contributes to the MLTA tendency. Among the zonal advection terms, the advection of mean temperature by anomalous current ( −u �T ∕ x ) is in phase with the MLTA tendency (Fig. 5c). This anomalous zonal advection term is − 0.04 °C month −1 , larger than that (− 0.004 °C month −1 ) in the Ctrl. Figure 6 shows the composite mixed layer terms during the period of the strong intensification of La Niña (Sep (2) to Dec (2)). The MLTA tendency in Vol is − 0.08 °C month −1 , significantly greater than that (− 0.03 °C month −1 ) in Ctrl (at the 95% confidence level on a Mann-Kendall test). The anomalous zonal advection term ( −u �T ∕ x ) and the thermocline feedback ( −w T � ∕ z ) and are still the main contributors to the decreasing trend of MLTA in Vol.
The contribution of the surface heat flux terms over the eastern Pacific is checked in Fig. 7a, c. During the period from Dec (1) to Mar (2), cooler surface heat flux occurs ( Hflx v < 0 < Hflx c , subscript 'v' and 'c' represent the results in Vol and Ctrl respectively) over the Niño 3.4 region (Fig. 7a), and shortwave radiation is the main control in Vol ( Sw ′ v < Sw ′ c < 0 ), causing the decreased SST. This is because the volcanic aerosols can directly reduce the surface shortwave radiation by scattering some solar radiation back into space at the top of atmosphere (Robock 2000). In the later part of year (2) (Fig. 7c), the reduced surface heat flux in Ctrl is slightly larger than that in Vol ( Hflx c < Hflx v < 0). The latent heat occupies the greatest contribution of the heat flux change in Vol

Role of the anomalous WNP anticyclone in determining the La Niña development
It has been suggested that the establishment of the anomalous Philippine Sea anticyclone plays an important role in the transition from an El Niño to a La Niña phase (Wang  1999,2000). Nevertheless, it is still unclear whether this anomalous anticyclone has an important influence on the warm-to-cold ENSO transition or not after large tropical volcanic eruptions; that is, do atmospheric processes play a role in the rapid decay of El Niño? To explore this, the correlation coefficient between the winter sea level pressure (SLP) anomalies (Dec (1) to Feb (2)) over the WNP and the Niño-3.4 SST anomalies during the La Niña mature phase (Dec (2) to Feb (3)) in Vol (25 cases) were calculated. The results show a significant correlation ( r = −0.56, p < 0.01 ) between them, meaning that the anomalous SLP over the WNP plays a critical role in the formation of La Niña. Figure 8 shows the evolution of composite SLP and 925 hPa wind anomalies during the warm-to-cold ENSO transition year. During the El Niño mature winter (Fig. 8a,  b), the WNPAC is clear in Vol, but not clear in Ctrl. The stronger easterly anomalies form on the south of the WNPAC and the cold Kelvin waves should be seen on the equator (Suarez and Schopf 1988;Battisti and Hirst 1989;Wang et al. 1999), rapidly decreasing the SST in Vol. The positions of the positive SLP anomaly center and the anticyclone anomalies in Vol are located further to the east than in Ctrl. This can enhance the equatorial easterly anomalies and accelerate the SST decreasing rate. In the following spring  (Fig. 8c, d), the anticyclone anomalies propagate slowly eastward, in agreement with Wang and Zhang (2002). In the autumn, the anticyclonic anomalies disappear, but the stronger anomalous easterlies occur over the equatorial Pacific in Vol (Fig. 8g, h).
The shortwave radiation decreases (increases) over the eastern (western) Pacific during the advanced decaying phase of El Niño in Vol (Fig. 7a, b). The eastward positioning of the WNPAC reduces precipitation and clouds, increasing the surface shortwave radiation ( Sw ′ v > Sw ′ c ). This can strengthen the east-west temperature gradient and gradually enhance the anomalous equatorial easterlies. The stronger equatorial easterly anomalies during the late part of year (2) in Vol (Fig. 8h, j) can enhance the total wind speed, increasing the upward latent heat flux (Fig. 7c) and cooling the equatorial eastern Pacific SST. Over the western Pacific, a largely warmer surface heat flux is found in Vol than in Ctrl, strengthening the anomalous equatorial easterlies.
To see the impact of the equatorial wind field change on ocean dynamics, the surface wind stress anomalies and the corresponding thermocline depth anomalies are shown in Fig. 9. The strong anticyclonic wind stress curl occurs in the equatorial western Pacific in winter (Fig. 9b), and moves eastward the following spring (Fig. 9d). The stronger equatorial easterly wind stress anomalies are enhanced on the south of WNPAC, uplifting the thermocline over the equatorial central and eastern Pacific in Vol. This can enhance the uplifting of thermocline and increase the anomalous vertical gradient ( , while the mean vertical upwelling term w v is almost the same as w c ( > 0 ). Hence, thermocline feedback is enhanced (Fig. 5b), leading to the larger downward MLTA tendency (Fig. 5a). Meanwhile, the anomalous westward movement of the mixed layer zonal current is strengthened ( u ′ v < u ′ c < 0 ), and the zonal mean temperature gradient term T v ∕ x is almost equal to T c ∕ x ( < 0 ), which can also contribute to the greater cooling rate (1) to Mar (2) (Fig. 5c).
During the development phase of La Niña (Fig. 9g-j), the WNPAC disappears, but the anomalous easterlies are stronger in Vol ( u ′ v < u ′ c < 0 ) which strengthens the anomalous zonal advection term . 6c) and cools the MLTA. Additionally, during the period of Sep (2) to Dec (2), the anomalous thermocline depth in Vol is still shallower than that in Ctrl ( T v � ∕ z > T � c ∕ z > 0 ), contributing to the greater cooling thermocline feedback . As the advanced decay of warming and the strong intensification of cooling over the equatorial eastern Pacific in Vol are related to the eastward positioning of the WNPAC, it is necessary to address how volcanic eruptions cause this.

Roles of the volcanic forcing in the eastward positioning of the WNPAC
The cause of the eastward shift of the WNPAC during DJF (1/2) in Vol (Fig. 3a, b) can be traced back to El Niño's development in year (1). During the developing phase of El Niño (Jun (1) to Nov (1)), strong warming occurred in the equatorial western-central Pacific (WCP, 5°S-5°N, 160°E-180°) in Ctrl, but little warming occurred in Vol (Fig. 10a, b). El Niño's warming also shifted eastward in Vol compared to Ctrl. To examine the process leading to El Niño's eastward shift in Vol, the MLTA budget analysis was used over WCP during the developing phase of El Niño (Jun (1) to Nov (1)) ( Table 2). The MLTA shows an increasing tendency in Ctrl ( T � c ∕ t = 0.12 °C month −1 ) during the developing phase of El Niño, and a slightly decreasing tendency in Vol ( T � v ∕ t = −0.01 °C month −1 ). The anomalous zonal advection term ( −u �T ∕ x > 0 ) brings heat to the mixed layer of WCP in both Vol and Ctrl. However, the differences in the dynamic term are primarily due to the vertical advection terms ( −w �T ∕ z and −w � T � ∕ z ), while the difference in the thermodynamic term is primarily attributed to the anomalous latent heat flux (Lh′).
The surface temperature anomalies and surface wind stress anomalies are shown in Fig. 11. Under tropical volcanic forcing, the cooling of land surface over the tropical area is greater than that over the ocean (Fig. 11a,  b). This can induce a stronger downward motion over the maritime continent, which pushes convective anomalies eastward (Fig. 11b). At the same time, the associated strong anomalous cyclonic wind stress curl over the CWP induces upward movement in the Ekman layer. This results in an eastward propagation of the anomalous oceanic upwelling to the CWP ( W v > 0 ) in Vol compared with that in Ctrl, causing the local cooling tendency of MLTA via the Ekman pumping term ( −w v �T v ∕ z < 0 ) over the CWP (Fig. 11d). After volcanic eruptions, the anomalous vertical temperature gradient term is larger , causing a more vigorous non-linear vertical advection ( −w � v T v � ∕ z < 0 ). The equatorial wind stress anomalies over the CWP can also enhance the total wind speed and the upward surface latent heat flux (Fig. 11b) . Therefore, the vertical advection terms and the upward surface latent heat flux suppress the warm SST over the CWP and then cause the eastward position of El Niño in the Vol.

Conclusions and discussion
Previous studies found that El Niño preferably occurs after tropical volcanic eruptions. This study presents evidence of the preferred La Niña occurrence in the second year after large tropical volcanic eruptions (year (2)) by examining the ensemble of nine sets of cold-season ENSO reconstructions. The processes and mechanisms of the warm-to-cold ENSO transition after tropical eruptions are investigated using CESM, offering the following major findings: 1. An analysis of the ensemble of nine sets of cold-season ENSO reconstructions reveals that a significant La Niñalike state takes place in year (2) after the large tropical volcanic eruptions. The numerical simulation shows that El Niño experiences an advanced decay and results in a stronger La Niña in Vol than the intrinsic ENSO mode in Ctrl during year (2). 2. The eastward positioning of the WNPAC in Vol plays a critical role in the earlier decay of warming and the strong intensification of cooling over the equatorial eastern Pacific (Fig. 12b). In these two periods, the easterly anomalies are enhanced on the southern edge of the anomalous anticyclone, which stimulates oceanic upwelling Kelvin waves, shallowing the thermocline and leading to the greater cooling rate by enhancing thermocline feedback and zonal advection feedback. Among the surface heat flux change terms over the equatorial eastern Pacific, the reduced shortwave radiation influences the earlier decay of warming, while the upward latent heat flux influences the strong intensification of the cooling.  3. The eastward positioning of volcanic-induced El Niño warming leads to the eastward positioning of the WNPAC (Fig. 12a). The eastward positioning of El Niño warming results from the volcanic aerosol effects. The volcanic aerosol effects cause stronger cooling over the maritime continent, causing an eastward shift in convective anomalies in Vol. At the same time, the strong anomalous cyclonic wind stress curl over the CWP enhances vertical thermal advection and upward surface latent heat flux, suppressing the development of warm SST over the CWP and causing the eastward positioning of the El Niño warming in year (1) in Vol.
Note that the volcanic forcing used in this study was reconstructed by Gao et al. (2008) who found that stronger cooling was produced in their model, such as the greater cooling which occurred in the year 1259, than in other reconstructions. The reason for this may be that the Ice-core Volcanic Index 2 overestimates the forcing of large volcanic eruptions owing to the linear assumption made, which is also possibly due to the uncertainties in both the model and reconstructions. This can cause cooling over the Niño-3.4 region in the immediate cold season after the large tropical eruptions (Fig. 1k). In addition, there are some uncertainties in different ENSO proxy data, although they correlate significantly with the modern observational data (Table 1). These uncertainties may cause some mismatches over the longer time period and the reason for this requires further exploration. Thus, the ensemble of nine cold-season ENSO reconstructions have been used in this study to reduce the uncertainties involved. Further, the SEA has been used to extract the composite ENSO and related atmospheric circulation anomalies after volcanic eruptions in Vol. However, it may not be possible to completely remove the internal variability during each volcanic eruption event because the initial conditions before the volcanic eruptions differ from each other. Therefore, multi-ensemble volcanic forcing experiments need to be conducted in future studies.

Fig. 12
Schematic diagram of the mechanisms responsible for La Niña's development after large volcanic eruptions. The elements in a, b are located at the sites where they take place. The blue shading denotes the reduced temperature, while the red shading represents the increased temperature. The gray shading represents the land distribution