Mechanisms linking multi-year La Niña with preceding extreme El Niño

El Niño/La Niña, characterized by anomalous sea surface temperature warming/cooling in the central-eastern equatorial Pacic, is a dominant interannual variability with irregularity, impacting worldwide weather and socioeconomics. The observed records show that La Niña often persists for more than two years, called “multi-year La Niña” which tends to accompany extreme El Niño in the preceding year; however, the physical linkage between them remains unclear. Here we show using reanalysis data that an extreme El Niño excites atmospheric conditions that favor the generation of the multi-year La Niña in subsequent years. Easterly wind anomalies along the northern off-equator in the Pacic during the decay phase of an extreme El Niño are crucial. They act to discharge ocean heat content (OHC) via an anomalous northward Ekman transport; the negative OHC anomaly is large enough to be restored by a single La Niña and, therefore, causes another La Niña to occur in the second year. Furthermore, analyses of the Coupled Model Intercomparison Project Phase 6 (CMIP6) models show that the occurrence frequencies of multi-year La Niña and extreme El Niño are highly correlated, supporting the abovementioned mechanism. Our results provide physical evidence that the increasing frequency of multi-year La Niña is explained by the increasing El Niño amplitude since the late 20th century.


Introduction
The El Niño-Southern Oscillation (ENSO) is a quasi-periodic variability with irregularity in several aspects such as amplitude and transition, and the life cycle is known to be locked with seasonal march: growing in boreal spring and summer, maturing in winter, and decaying in the following spring [1][2][3] . The oscillatory nature of ENSO has been well explained by the linear theory called the recharge oscillator paradigm, in which equatorial Paci c warm water volume anomalies can be used as a harbinger of the ENSO sea surface temperature (SST) anomaly [4][5][6] . According to this theory, the transition from La Niña event to El Niño event occurs such that the intensi ed easterly trade winds associated with the La Niña SST anomaly act to recharge warm subsurface water from the off-equator to the equatorial strip via a slow Sverdrup adjustment. The deepened thermocline weakens the cooling effect due to mean upwelling, referred to as the thermocline feedback, which diminishes the La Niña SST anomaly and, then, causes ENSO turnabout. This process has been veri ed using a phase relationship between the observed warm water volume anomaly in the equatorial Paci c and the Niño 3.4 SST anomaly [7][8][9] .
In reality, the temporal behavior of ENSO is highly complex 10 . Unlike solutions in the recharge oscillator equations, observations show that ENSO has a strong transition asymmetry; El Niño terminates rapidly but La Niña often lasts more than two years, the latter referred to as multi-year La Niña 11 . In particular, extremely strong El Niño events tend to terminate and transition to La Niña within one year. When extreme El Niño decays from winter to spring, equatorial westerly wind anomalies are known to shift southward 12 , causing equatorial thermocline shoaling and contributing to the rapid termination of extreme El Niño [13][14][15] . This atmospheric pattern with southward-shifted westerly anomalies is induced by a nonlinear coupling between ENSO's interannual frequency and the annual cycle in the western Paci c, known as a combination mode (C-mode) [16][17][18] .
The multi-year persistency of La Niña has been investigated using observations and climate model simulations. A key process for the ENSO phase transition is the delayed negative thermocline feedback (equivalent to the recharge/discharge process), which is weak for La Niña and may, therefore, explain the longer duration of La Niña than El Niño 19 . A report also states that the SST anomaly pattern has a wide meridional structure when La Niña persists for more than two years 20 . Because the recharge after the peak of La Niña is effective with a narrow meridional structure of the wind stress curl, the difference in the SST pattern may contributes the multi-year persistency of La Niña. In addition, remote in uence from the Atlantic and Indian Ocean may play a role in the multi-year persistence of La Niña that was predicted out to two years ahead in a model 21 . A common suggestion from previous studies is that a multi-year La Niña tends to occur after an extreme El Niño. A recent study shows that the duration of La Niña is strongly in uenced by the amplitude of the preceding El Niño in both observations and a long climate model simulation, presumably due to a large initial discharge 22 . However, the physical link between extreme El Niño and the subsequent multi-year persistence of La Niña has not been comprehensively examined using ocean subsurface reanalysis data. Furthermore, extremely strong El Niño events decay rapidly owing to a coupling with the annual cycle 17 ; however, the in uence of the rapid termination of El Niño on the occurrence of multi-year La Niña is unclear.
In this study, we address the following questions: What process distinguishes multi-year and single-year La Niña? Is there a dynamical connection between extreme El Niño and a subsequent multi-year La Niña? By answering these questions, we identi ed triggers for multi-year La Niña. In addition to the observational analyses, we investigated the connection between extreme El Niño and multi-year La Niña in 37 CMIP6 models and demonstrated that their observed physical linkage is consistently seen in the climate model ensemble. There are several reports that a multi-year La Niña in the rst and second years has distinct impacts on the atmospheric circulation over North America and East Asia 11,23,24 . Although ENSO forecasting for more than one year is still challenging 25,26 , deepening our understanding of the mechanism of a multi-year La Niña and its connection to a preceding extreme El Niño would help improve seasonal and long-term climate prediction.

Time evolution of two types of La Niña
We identify six and four multi-year and single-year La Niña, respectively, using the observed Niño 3.4 SST anomaly (hereafter N3.4) during 1961-2016 ( Supplementary Fig. 1). We de ne the period when the rst La Niña develops as Year (0) in the composite analysis. By taking composites for these events of N3.4 and the ocean heat content (OHC) integrated over the equatorial Paci c (OHC eq : 120° E-60° W, 5.5° S-5.5° N, surface-500 m), the latter measuring the equatorial Paci c warm water volume, trajectories associated with the two types of La Niña are obtained in the phase space (Fig. 1a). The composite-mean trajectories are akin to the limit cycle 4,7 , showing that La Niña in Year (0) follows El Niño in Year (-1), both of which reach their peak in boreal winter (November-December-January: NDJ, circles in Fig. 1a). Four out of six multi-year La Niña events accompany extreme El Niños (the NDJ-mean N3.4 exceeding the 1.5 standard deviation) in the preceding year, whereas all El Niño before a single-year La Niña has a moderate amplitude ( Supplementary Fig. 1). Consistent with the recharge theory, the OHC eq anomaly lead N3.4 by 2-3 months and peaks in September-October-November (SON) before the mature phase of La Niña. While the positive OHC eq anomaly during SON(-1) is comparable in magnitude between the two composites, the negative OHC eq anomaly in SON(0) preceding multi-year La Niña is four times as large as that of single-year La Niña. Despite the large difference in the OHC eq anomaly, the negative N3.4 in NDJ(0/1) in the multi-year composite is only 1.3 times larger than that in the single-year composite due to a strong asymmetry in the thermocline feedback and nonlinear dynamical heating, which suppresses La Niña but enhances El Niño 7,27,28 . This contrast in magnitude between the SST and OHC eq anomalies is also an important feature of multi-year La Niña and may contribute to the recurrence of La Niña.
After NDJ(0/1), the negative OHC eq anomaly eventually turns positive in the single-year composite, leading to the termination of La Niña. In contrast, the OHC eq anomaly in the multi-year composite decays but remains negative until the following year, resulting in an intensi ed subsequent La Niña in NDJ(1/2).
The time evolution shown in Fig. 1a indicates that the recurrence of La Niña is consistent with the linear recharge oscillator theory. An exceedingly negative OHC eq anomaly in NDJ(0/1) that is not restored within a year is crucial in understanding multi-year La Niña 22 . The spatial distribution of composite OHC anomalies in its peak season, SON(0), is shown in Fig. 1b,c. The multi-year composite is characterized by a contrast between large negative OHC anomalies in the equatorial Paci c and positive anomalies in the western-central Paci c at approximately 5°-15° N, suggesting that the former was induced by a mass exchange between the equatorial strip and northern off-equator. Such an equatorial asymmetry is not observed in the single-year composite, which indicates weak negative OHC anomalies in the centraleastern equatorial Paci c and positive anomalies in the tropical western Paci c.
Physical link between multi-year La Niña and preceding extreme El Niño A signi cant difference between multi-year and single-year La Niña is clearly observed in the composite time series of N3.4 and OHC eq anomalies ( Fig. 2a,b). To clarify the processes responsible for the strong discharge of warm water during the growth phase of multi-year La Niña, heat budget analysis was performed on OHC eq using ocean reanalysis datasets (Methods). We present the result based on ORAS4 in Fig. 2b-d because the budget terms reveal the smallest error among the four reanalysis data; however, our conclusion depends little on the dataset ( Supplementary Fig. 2).
A large negative OHC eq anomaly in SON(0) for the multi-year La Niña composite occurs because of a large negative OHC eq tendency during the rst half of Year (0), i.e., from February to June, denoted as FMAMJ(0) (black solid curve in Fig. 2c). For both multi-year and single-year La Niña composites, the OHC eq tendency is reproduced well by the aggregated OHC eq budgets (shading in Supplementary  Fig. 2a,e), which can then be decomposed into individual terms (zonal and meridional advection terms due to geostrophic and Ekman currents, heat exchange at the bottom of the subsurface, and net heat ux at the surface). The sum of the geostrophic and Ekman terms (i.e., Sverdrup heat transport) is called the recharge rate, which explains most of the OHC eq tendency (green curves in Fig. 2c). Unlike the linear recharge oscillator theory, the composited recharge rate indicates a persistent large negative anomaly during FMAMJ(0) when the N3.4 signal is weak.
The decomposition of the recharge rate into components of surface Ekman current heat transport (EHT) and geostrophic current heat transport (GHT) is shown in Fig. 2d. On the one hand, the GHT during the peak period of the preceding El Niño shows a large negative anomaly in the multi-year composite compared to the single-year composite; however, they are similar after scaling with N3.4 (Supplementary Table 1). This implies that the GHT-induced discharge is not a primary factor for generating the OHC eq difference between multi-year and single-year La Niña. On the other hand, EHT in the multi-year composite becomes negative in FMAMJ(0) and contributes to the large negative recharge rate (i.e., discharge); hence, a large negative OHC eq tendency is not observed in the single-year composite (Fig. 2c).
The composite heat transport at the northern and southern boundaries indicates that the strong discharge for multi-year La Niña in FMAMJ(0) is attributed to EHT across the northern boundary, which is slightly weakened by the heat imported across the southern boundary ( Supplementary Fig. 3). In general, GHT is inversely proportional to N3.4 throughout the year and primarily represents the recharge theory, while EHT has a similar magnitude to GHT but is proportional to N3.4 from the summer to late autumn, indicating that this process disturbs a theoretical recharge/discharge ( Supplementary Fig. 4c-f). Consequently, the recharge rate is not correlated with N3.4 during the summer (Supplementary Fig. 4a,b).
During FMAMJ(0) corresponding to the transition period from El Niño to La Niña, composite anomalies of SST, precipitation, and surface wind stress show different patterns between multi-year and single-year La Niña (Fig. 3). The El Niño SST pattern persists in the multi-year La Niña composite but transforms to a weak La Niña pattern in the single-year composite (Fig. 3a,c) because strong El Niño tend to last several more months compared to weak El Niño 22 . In the multi-year composite map, anomalous positive precipitation is found to the south of the equator, while negative precipitation anomalies appear along the northern off-equator. This southward shift of the precipitation anomaly pattern has been identi ed during the decay phase of extreme El Niño and is because the climatological mean SST distribution shifted southward in this season 16 . Consistently, eastward wind stress anomalies are present over the south of the equator, whereas westward anomalies dominate over the northern off-equator, the latter being responsible for the strong EHT-induced discharge (Fig. 3b). The anomalous cross-equatorial wind stresses can form in response to diabatic heating anomalies associated with the southward-shifted precipitation anomaly pattern shown in Fig. 3a (Supplementary Fig. 5a,b). Overall, anomaly patterns in precipitation and wind stress in the multi-year composite are consistent with the equatorial asymmetric distribution of OHC anomalies during SON(0) (Fig. 1b). The atmospheric responses are similar to the Cmode, which acts to terminate extreme El Niño 13,18 , supporting links with a preceding extreme El Niño.
The single-year composite maps lack anomalous zonal wind stresses over the off-equator, likely due to weak, symmetric precipitation anomalies around the equator (Fig. 3c, d).
Exceptional cases of coupling between extreme El Niño and multi-year La Niña exist. Yet, the importance of the easterly wind stress anomaly over the northern off-equator for multi-year La Niña is commonly identi ed. Two out of six multi-year La Niña events (1970-72 and 2007-09) are preceded by moderate but not extreme El Niño (Supplementary Fig. 6). Composite SST and precipitation anomalies for these two events are almost similar to those of the composites of single-year La Niña. However, the anomalous easterly wind stresses appear over the northern off-equator as a response to thermal forcing in the equatorial Paci c (Supplementary Fig. 5c,d; Supplementary Fig. 6c, d). In addition, two out of six extreme El Niño events (1965-66 and 1991-92) do not accompany a multi-year La Niña ( Supplementary Fig. 7).
The composite anomalies for these two El Niño events show that the negative OHC eq anomaly is weak in Year (0/1). Consistently, easterly wind stress anomalies are absent along the northern off-equator during FMAMJ(0). These results con rm that the anomalous easterly winds and associated EHT-induced discharge in the northern off-equatorial latitudes link multi-year La Niña with a preceding El Niño.

Analysis of CMIP6 control simulations
To obtain a robust relationship between extreme El Niño and multi-year La Niña, we repeated our analyses in 500-year long preindustrial control simulations using 37 Earth system models (ESMs) participating in CMIP6 29 (Supplementary Table 2 year La Niña occur less frequently in ESMs than in observations; however, some models capture the observed frequency of multi-year events. As expected, a signi cant positive correlation ( ) is found between the frequencies of multi-year La Niña and extreme El Niño in CMIP6 models (Fig. 4a). This indicates that multi-year La Niña tend to occur more often in a model that generates many extreme El Niño, supporting a physical link between them.
We calculate the ratio of the number of multi-year La Niña events accompanying preceding extreme El Niño against the total number of multi-year La Niña events. This ratio, denoted as R E−L , measures the coupling between extreme El Niño and multi-year La Niña (Fig. 4b) for the entire ensemble). Indeed, the average ENSO amplitude in ten highest R E−L models is approximately 50% larger than the average amplitude in ten lowest R E−L models ( K versus K). This relationship can be understood given that a frequent occurrence of extreme El Niño which increase both the overall ENSO amplitude and frequency of multi- year La Niña. Thus, the analysis of CMIP6 models supports the physical mechanism linking extreme El Niño with multi-year La Niña. An increase in observed R E−L in the late 20th century may also be related to an increase in ENSO amplitude.

Implications for ENSO dynamics
In this study, we examined the mechanisms of multi-year La Niña which occupies two third of La Niña events during the period 1961-2016. A physical link was determined between multi-year La Niña and extreme El Niño, of which the latter often accompanied in the preceding year. The essential mechanism for the transition from an extreme El Niño to multi-year La Niña is explained by a modi cation of the recharge/discharge cycle, which is schematically illustrated in Supplementary Fig. 8. Using composite analyses of observations and reanalysis data, we found that signi cant negative OHC anomalies form in the equatorial Paci c during the rst La Niña and they persist after its peak. This strong mass discharge of the upper ocean cannot be restored by a single La Niña and, therefore, causes another La Niña to occur in the second year. The large negative OHC anomaly before the peak of the rst La Niña is induced by intensi ed northward Ekman heat transport driven by surface easterly anomalies along the northern offequator. The patterns of surface wind stress anomalies, easterly anomalies to the north, and westerly anomalies to the south of the equator have occurred when extreme El Niño decays. This plays a central role in linking extreme El Niño with subsequent multi-year La Niña.
Given that ENSO is positively skewed 27,28,30 , i.e., El Niño is stronger than La Niña, it is reasonable to doubt that multi-year La Niña is an apparent feature arising from the climatological mean shifted toward El Niño due to amplitude asymmetry. However, our results clearly show that multi-year La Niña is not a statistical artefact but a part of the intrinsic complex nature of ENSO. Previous studies have shown the contribution of anomalous wind stress curl causing GHT in the ENSO transition 9,20,31−33 , an active role of inter-basin interaction for the transition asymmetry 34 , and a meridional shift of westerly surface wind anomalies for the transition from extreme El Niño to La Niña [12][13][14]35 . Although we do not exclude these processes during ENSO phase transition, we demonstrate that the anomalous Ekman heat transport is crucial to transition from extreme El Niño to multi-year La Niña. The link between extreme El Niño and multi-year La Niña as represented by R E−L is robust; therefore, the occurrence of multi-year La Niña may be predictable beyond the typical predictability of ENSO events 36 . If a coupled atmosphere-ocean model initialized with an extreme El Niño condition can reproduce the wind stress pattern responsible for the anomalous Ekman heat transport, the subsequent two years, when a multi-year La Niña will occur, maybe predicted owing to the large memory in the ocean heat content 21,37

Observations and reanalysis data
We used the observed monthly SST and precipitation datasets from COBE-SST2 40 and PRECipitation REConstruction (PREC) 41 of the National Oceanic and Atmospheric Administration, respectively. Four ocean reanalysis products that contain temperature, salinity, current velocity, and surface heat ux were also used: Ocean Reanalysis System 4 (ORAS4) 42 , Ocean Reanalysis System 5 and its backward extension (ORAS5) 43 , German contribution to the consortium for Estimating the Circulation and Climate of the Ocean system 2 (GECCO2) 44 Fig. 1). All La Niña events were preceded by El Niño. When the multi-year La Niña lasted for three years, we analyzed the rst two years. For the CMIP6 multi-model analysis, the threshold to detect El Niño/La Niña was 0.7 standard deviation, which is 0.5 K in COBE-SST2 for the period 1901-2018, but corresponds to a different value in a different model.

Diagnosing steady atmospheric response to an imposed heating
We diagnosed a thermally forced steady atmospheric response using a linear baroclinic model (LBM) 49 . The model has a horizontal resolution of T42 and 20 vertical levels. The basic states, FMAMJ climatology, and anomalous diabatic heating (including condensational heating, radiative heating, surface heat ux, and subgrid-scale heat ux convergences) 50 were calculated using the combined data of ERA40 and ERA-Interim. For simpli cation, the thermal forcing was con ned to the tropical Paci c (120° E-80° W, 30° S-30° N).

Statistical signi cance
We applied a two-sided Student's t-test for the composite and correlation analyses, and Welch's t-test for the difference between multi-year and single-year composite means. The con dence levels are described in the gure captions.