Rainfall sustains multiyear La Niña

Domain: earth_science
Published: 2026-01-14
Editorial outcome: minor_revision

Abstract
Historical observations and climate models indicate that multiyear La Niña events are growing more frequent. These prolonged events cumulatively alter tropical Pacific rainfall, which in turn affects regional mixed-layer salinity (MLS). However, how cumulative rainfall anomalies modulate multiyear La Niña through salinity effects remains unclear. Here, using observations and model experiments, we show that the positive MLS anomalies in the western–central equatorial Pacific—driven initially by oceanic dynamical processes and later by rainfall deficits during multiyear La Niña—reinforce these prolonged La Niña events. Quantitatively, the rainfall reduction enhances La Niña amplitude by 14% in the first year and 32% in the second year. MLS anomalies in the western–central equatorial Pacific initially trigger equatorial Kelvin wave adjustments, causing rapid surface cooling in the eastern equatorial Pacific. Subsequent slow ocean circulation responses lead to basin-wide equatorial Pacific cooling. The superposition of rapid and slow oceanic responses creates cumulative positive feedback from rainfall on the second-year La Niña. These findings identify rainfall–salinity feedbacks as a key mechanism sustaining multiyear La Niña events, with implications for ENSO prediction and hydrological cycle.

Introduction
Our understanding of the diversity of the El Niño-Southern Oscillation (ENSO) has improved considerably over the past four decades1–4, as evidenced by the various types of El Niño events, the amplitude asymmetry between El Niño and La Niña, and the tendency for La Niña to persist for more than 2 years, known as multiyear La Niña events5–9. Recent studies have indicated a marked increase in the occurrence frequency of multiyear La Niña events during the historical period10, and climate models project a further significant rise in their frequency throughout the 21st century11. Consequently, the importance of multiyear La Niña events to global weather and climate is garnering increased attention from the scientific community7,8,12.

Recent studies on the formation and development of multiyear La Niña events have predominantly focused on ocean–atmosphere interactions13–18. A widely accepted hypothesis is that multiyear La Niña events are often preceded by a strong El Niño19. The surface wind anomalies associated with a strong El Niño can lead to a significant discharge of heat content in the equatorial ocean, requiring the tropical Pacific coupled system to take more than a year to return to a neutral state. Studies have also emphasized the impacts of the mid-latitude Pacific Ocean and other basins on multiyear La Niña events, particularly the North Pacific20,21. In addition, the persistence of multiyear La Niña events can be attributed to a delayed heat content recharge driven by anomalous Ekman transport, nonlinear thermocline feedbacks, and interactions with the Indian and Atlantic Oceans13,22–24. Nevertheless, during multiyear La Niña events, the equatorial ocean undergoes sustained cooling for over 2 years—a timescale that implies slow ocean dynamical processes (e.g., vertical mixing, wind-driven large-scale circulation, and rainfall-induced salinity changes) may play a pivotal role.

Among these slow ocean processes, salinity variations induced by rainfall play a critical and often overlooked role. Rainfall fundamentally modulates ocean salinity in the tropical Pacific across a range of timescales, from the mean state to interannual variability25–30. During ENSO events, the mixed layer salinity (MLS) exhibits a well-defined spatial structure that closely corresponds to anomalous rainfall patterns30. This interannual variability is further influenced by zonal advection of low-salinity waters, driven by anomalous eastward surface currents during El Niño events. Conversely, during La Niña, anomalous westward currents lead to the opposite effect31–34. Additionally, rainfall variability—tied to shifts in atmospheric deep convection and the zonal migration of the western Pacific warm pool—exerts a significant influence on these salinity patterns34–36.

In turn, MLS variability associated with rainfall fluctuations during the ENSO cycle profoundly affects upper-ocean stratification, mixed layer dynamics and hence temperature33,37. For example, accounting for salinity effects significantly improves the simulation of the equatorial Pacific cold tongues38. Including salinity effects led to a 0.8°C increase in sea surface temperature (SST) near the eastern boundary of the warm fresh pool during the 1997–1998 El Niño36. The salinity-induced barrier layer also plays a critical role in El Niño development by suppressing the upward entrainment of cold subsurface water into the mixed layer39–41. Rainfall-induced salinity anomalies also have an important impact on ENSO diversity. For instance, rainfall introduces a positive feedback mechanism that amplifies ENSO variability42. Through this feedback, differing zonal patterns of salinity anomalies enhance ENSO asymmetry (manifested as larger SST anomalies during El Niño) and also account for ~10% of the amplitude difference between Eastern Pacific (EP) and Central Pacific (CP) El Niño events43,44.

Therefore, it is hypothesized that since rainfall-induced MLS anomalies can persist for over 2 years during multiyear La Niña events, they may exert a sustained influence on the upper ocean, which, in turn, could sustain the cooling. The underlying processes are plausible: the absence of freshwater input from negative anomalies in rainfall reduces upper-ocean stratification, thereby enhancing vertical mixing and entrainment, deepening the mixed layer, and modulating large-scale ocean circulation and equatorial waves. These processes collectively regulate SST evolution over extended timescales in a manner that favors the persistence of La Niña. However, the extent to which these rainfall-driven MLS variations govern the progression of multiyear La Niña events remains undetermined. In this study, we address this issue using observational data and numerical simulations, demonstrating that multiyear La Niña events exhibit a self-amplifying nature through a cumulative rainfall–salinity feedback mechanism.

Results
Asynchronous MLS responses to rainfall and ocean dynamics
Interannual anomalies of rainfall and MLS are closely tied to ENSO evolution. In the Niño4 region (5°S–5°N, 160°E–150°W), interannual variations in MLS and rainfall exhibit strongly synchronized fluctuations, with a Pearson correlation coefficient of −0.68 (statistically significant at the 95% confidence level, two-tailed t-test). This robust negative correlation underscores a strong linkage between rainfall deficits and positive MLS anomalies in the western-central Pacific. To characterize the MLS-rainfall relationship during multiyear La Niña events, we selected six such events between 1980 and 2022 for composite analysis (Fig. 1a). Two of these events (1998–2000 and 2020–2022) persisted into a third year while maintaining a cooling state. Throughout these events, persistent negative rainfall anomalies coincided with positive MLS anomalies in the western-central Pacific, though their amplitudes varied—likely reflecting differences in event intensity. A pronounced negative rainfall anomaly was observed near the dateline, associated with persistent surface cooling in the western-central equatorial Pacific (Fig. 1b, c). Climatologically, low-salinity waters (e.g., <34.8 psu) dominate the western-central equatorial Pacific, sustained by abundant rainfall. During multiyear La Niña, the isohaline shifts westward by ~10°, indicating a westward expansion of higher-salinity water (Supplementary Fig. 1). Together with concurrent rainfall deficits, this zonal migration of isohalines suggests that zonal advection contributes to MLS variability during certain phases of multiyear La Niña. Similar features are reproduced across multiple reanalysis and observational datasets (Supplementary Fig. 1).

Nevertheless, the specific roles of rainfall versus oceanic dynamical processes during these events remain unclear. To address this, we employed an MLS budget analysis using ORAS5 data to quantify the contributions of different processes. We focused on three key periods for the Niño4 region: MAMJ(0), ND(0)JF(1), and ASON(1). These represent, respectively, a rapid MLS increase during the first-year La Niña, the MLS peak within that year, and a second MLS peak late in the second-year La Niña (Supplementary Fig. 2). Here, “(0)” denotes the first year and “(1)” the second year of a multiyear La Niña event. The MLS budget analysis in the Niño4 region reveals that the dominant processes contributing to salinity anomalies vary by phase: (a) During MAMJ(0), as the first-year La Niña develops, surface waters in the Niño4 region undergo salinification primarily driven by ocean dynamical processes, notably zonal advection and vertical mixing associated with mixed layer deepening (Fig. 2a). In contrast, the contribution from freshwater flux (E-P) is minor. (b) During ND(0)JF(1), at the first MLS peak, a distinct rainfall reduction emerges. This reduction in rainfall results in a significant positive contribution from freshwater flux to MLS anomalies. Conversely, horizontal advection counteracts this increase, driving a temporary decrease in MLS. (c) During ASON(1), at the second MLS peak, positive MLS anomalies are again primarily attributable to reduced rainfall.

Zonal advection anomalies contribute robustly to increasing MLS in the Niño4 region during the first-year La Niña. This is primarily due to the advection of mean salinity (S̅) by anomalous zonal currents (u′), represented by the term −u′∂S̅/∂x (Fig. 2b–e). Westward current anomalies (u′ < 0)—a wind-driven response to La Niña–related easterly wind anomalies (Supplementary Fig. 3)—act on the mean zonal salinity gradient (∂S̅/∂x > 0), promoting salinification in the Niño4 region via the term −u′∂S̅/∂x, which is positive during MAMJ(0). A comparable mechanism operates during subsequent peak phases (Fig. 2b; ND(0)JF(1) and ASON(1)). The contribution of meridional advection to MLS variations is relatively small and is linked to the migration of the intertropical convergence zone. In summary, reduced rainfall significantly contributes to MLS increases during the peak of the first-year La Niña and continues to modulate MLS evolution in the second year, whereas zonal advection and vertical mixing play a pivotal role in the initial MLS increase during the first-year La Niña. Thus, the MLS responses to rainfall and ocean dynamic processes are asynchronous throughout the multiyear La Niña period.

Delayed MLS response to rainfall
To further investigate the effects of rainfall anomalies on MLS during multiyear La Niña events, we performed two sets of experiments using two different ocean circulation models: the Gent-Cane (GC) ocean general circulation model (OGCM)45 and the Regional Ocean Modeling System (ROMS)46. In the rain_inter experiment, rainfall was specified to vary interannually, whereas in the rain_clim experiment, it was prescribed as a climatological field; the difference between rain_inter and rain_clim represents the oceanic response to interannually varying rainfall.

The GC-OGCM effectively simulates the characteristic MLS response to multiyear La Niña events (Fig. 3a), including its delayed response to rainfall. In the rain_inter experiment, the onset of positive MLS anomalies during MAMJ(0) is closely aligned with that observed in the IAP and ORAS5 datasets (Supplementary Fig. 1). However, due to the weaker vertical mixing and entrainment associated with wind stirring in the GC-OGCM47, the simulated amplitude of interannual MLS anomalies is subsequently underestimated. When interannual rainfall anomalies are removed in the rain_clim experiment, the initial positive MLS anomalies are still effectively preserved (Fig. 3b). This indicates that dynamic processes (e.g., advection and vertical mixing) are the primary drivers of positive MLS anomalies during the developing phase of the first-year La Niña, which are subsequently sustained through the following autumn and winter. Figure 3c shows the difference in MLS between the rain_inter and rain_clim experiments, highlighting the direct impact of rainfall. The results demonstrate that the influence of interannual rainfall variability on MLS primarily manifests during the peak phase of the first-year La Niña and persists into the third year of the event, underscoring a delayed MLS response to rainfall forcing.

Similar results are reproduced in the ROMS simulations (Fig. 3d). The spatial evolution of the positive MLS anomaly aligns well with observational data, although its onset (around July–August of the first-year La Niña) is slightly delayed compared to both observations and the GC-OGCM simulations. The rain_clim experiment reveals that the MLS anomaly exhibits two prominent peaks during the autumns of the first- and second-year La Niña (Fig. 3e). Notably, positive MLS anomalies emerge earlier in the first year in the rain_inter than in the rain_clim, indicating that oceanic dynamic processes alone can initiate these anomalies. Regarding the effects of interannual rainfall, the ROMS results are consistent with those from the GC-OGCM, showing a clear delayed effect on ocean salinity (Fig. 3f). Moreover, rainfall continues to dominate the positive MLS anomalies from the first year into the second year of the multiyear La Niña (Fig. 3f).

In summary, results from both models highlight two key findings: oceanic dynamic processes are responsible for the initial development of positive MLS anomalies at the onset of multiyear La Niña, whereas rainfall plays a predominant role in sustaining and amplifying these anomalies during the subsequent evolution.

Effects of rainfall on ocean stratification during multiyear La Niña
Given the significant impact of rainfall on MLS during multiyear La Niña events, we next investigate how these rainfall-induced MLS anomalies, in turn, affect SST by modulating the upper-ocean structure. We analyze the influence of MLS changes on upper-ocean stratification and vertical mixing, and further explore the salinity-mediated pathway through which rainfall affects SST, using GC-OGCM experiments. Figure 4b, c show composite potential temperature and salinity anomalies during multiyear La Niña events, revealing negative temperature and positive salinity anomalies in the western–central equatorial Pacific. These anomalies collectively result in positive density anomalies beginning in the summer of the first-year La Niña (Fig. 4a).

The composite shows that negative anomalies in the squared buoyancy frequency (N²) develop below the mixed layer starting in September(0) and persist until the end of the event (Fig. 4d). This indicates a weakening of stratification in the western Pacific upper ocean due to a reduced vertical density gradient, which would enhance vertical mixing and entrainment at the mixed layer base41. Positive salinity anomalies play a crucial role in enhancing instability near the mixed layer base. When the effect of salinity is excluded, the composite N² anomaly shows a pattern similar to the full field, indicating that temperature is the primary driver of N² evolution (Fig. 4e). Nevertheless, the influence of salinity on N² becomes substantial near the base of the mixed layer. This occurs because salinity-induced density anomalies are concentrated within the mixed layer (Fig. 4c), thereby reducing the vertical density gradient (Fig. 4f). Given the key role of salinity in determining density at this depth, the resultant reduction in stratification enhances vertical mixing and entrainment, ultimately contributing to surface cooling.

Effects of rainfall on sustaining multiyear La Niña
To further investigate the cumulative effects of rainfall in sustaining multiyear La Niña, we conducted a sensitivity experiment using the GC-OGCM, in which the composite negative rainfall anomalies during these events were doubled. This adjustment was designed to compensate for the underestimation of MLS anomalies in the model (Fig. 3a) and to better isolate the salinity effect. The experiment yielded positive MLS anomalies that more closely matched observations in amplitude (Supplementary Fig. 4). Doubling the rainfall anomalies enhanced surface salinification and contributed to significant basin-wide cooling during multiyear La Niña events (Fig. 5a, b). The spatial patterns of these anomalies closely resembled those observed, particularly from the end of the first year through the second year (Supplementary Fig. 4). Once established, these positive salinity anomalies weaken the upper-ocean stratification, facilitating the entrainment of cold subsurface water into the mixed layer and leading to cold SST anomalies. This subsurface cooling peaked in the autumn and winter of the second-year La Niña across the equatorial Pacific (Fig. 5b). Quantitatively, the enhanced rainfall reduction amplified the La Niña amplitude by 14% in the first year and 32% in the second year, as measured by the Niño3.4 index during the respective winter peaks (Supplementary Table 1). In magnitude, despite the simulated MLS anomalies remaining weaker than observed, they still induced cold SST anomalies exceeding 0.3°C in the central-eastern equatorial Pacific.

Most studies have confirmed that imposing negative rainfall anomalies over the western tropical Pacific, with an amplitude comparable to observed interannual variability, can lead to a cooling exceeding –0.3°C in the central-western equatorial Pacific. This magnitude of cooling is consistent with the SST changes simulated in our experiment (Supplementary Table 2).

The specific role of rainfall-induced salinity changes in multiyear La Niña events has not been fully isolated. To better understand these dynamics, we performed an idealized experiment that mimicked the observed salinity anomalies. In this experiment, we prescribed a negative rainfall anomaly of −9 mm day⁻¹ within the region (140°E–180°, 5°S–5°N) from June(0) to December(1), designed to influence salinity alone (see the Experiments II in the Methods section for details). The resulting positive salinity anomalies were largely confined to the observed region (Fig. 1c) but were slightly stronger in magnitude (Fig. 5c and Supplementary Fig. 4). This idealized experiment demonstrates that the impact of rainfall-induced salinity on SST is consistent with the results from the doubled rainfall anomaly experiment (Fig. 5b). Positive MLS anomalies in the western-central equatorial Pacific induce local SST cooling. Although the MLS anomalies remain relatively constant, the resultant SST cooling gradually intensifies and spreads across the entire equatorial Pacific during the summer and autumn of the first-year La Niña. This effect leads to a sustained cooling into the second year (Fig. 5d).

We also conducted separate experiments applying the negative rainfall anomaly to the first-year (1yr_exp) and second-year (2yr_exp) La Niña conditions. The results indicate that the rainfall anomaly applied in the 1yr_exp has a greater impact on the overall cooling in the second year than that in the 2yr_exp, suggesting a cumulative effect of positive MLS anomalies on SST (Supplementary Fig. 5). Furthermore, additional experiments with rainfall anomalies applied at different longitudes confirm that the amplitudes of the SST and MLS responses during multiyear La Niña events are generally consistent across different spatial distributions of rainfall (Supplementary Fig. 6). These experiments collectively demonstrate that the effects of negative rainfall anomalies from the first and subsequent years superimpose, contributing to the more pronounced cooling in the second-year La Niña.

Mechanism of rainfall in sustaining multiyear La Niña
A mixed-layer heat budget analysis was performed to identify the controlling mechanisms in the Niño4 and Niño3 regions. In the Niño4 region (Fig. 5e), the imposition of negative rainfall anomalies in the idealized experiment induces a negative temperature tendency from June(0) to October(1), indicating a persistent cooling trend. This cooling is primarily driven by enhanced vertical mixing resulting from the rainfall-induced weakening of stratification at the base of the mixed layer (as shown in Fig. 4d). However, the cooling effect of vertical mixing is partially counteracted by the surface net heat flux. During the second-year La Niña, the contribution of zonal advection shifts from positive to negative, further enhancing the cooling in Niño4 region. In the Niño3 region, the mixed-layer temperature shows little change during the first-year La Niña due to the absence of significant positive salinity anomalies. In the second year, however, intensified vertical mixing causes the temperature tendency to become negative, a state that persists until the end of the year. Throughout this period, zonal advection, meridional advection, and heat flux damping alternately offset the cooling induced by vertical mixing (Fig. 5f). Therefore, two consecutive years of negative rainfall anomalies promote sustained equatorial Pacific cooling by modulating ocean salinity.

The oceanic dynamical response to idealized rainfall forcing occurs in two distinct stages: (1) a fast adjustment during the initial phase (June(0)–July(0)), and (2) a slow response driven by the cumulative effects of rainfall as the multiyear La Niña develops (Fig. 6a, b). In the fast response stage, the prescribed negative rainfall anomaly over the western-central equatorial Pacific increases mixed-layer salinity and density. This elevated density weakens upper-ocean buoyancy, leading to a deeper mixed layer. The deepening mixed layer triggers a high-order baroclinic adjustment, generating an alternating cold-warm-cold-warm temperature structure in the vertical dimension (Fig. 6c), whose vertical wavelength corresponds to the mixed layer depth48. These subsurface thermal anomalies propagate eastward as Kelvin waves; notably, a cold anomaly signal concentrated within the core of the Equatorial Undercurrent (EUC) propagates rapidly to the eastern equatorial Pacific via equatorial wave dynamics. This process leads to the emergence of a surface cooling signal in the eastern basin within 2–3 months, consistent with the propagation speed of Kelvin waves (Fig. 6d, e).

The slow response of oceanic processes stems from the persistence of positive MLS anomalies in the western-central equatorial Pacific. These anomalies maintain a deeper and denser mixed layer, which sustains enhanced vertical mixing and a cold anomaly in the region for up to two years (Fig. 6a, f). Furthermore, the positive MLS anomalies affect the sea surface height in the western-central equatorial Pacific. The resultant reduced zonal sea surface height gradient strengthens the westward surface currents; under intensified trade winds, this strengthens meridional advection and upwelling in the eastern equatorial Pacific, thereby reinforcing the basin-wide cold anomaly (Fig. 6f). Concurrently, an anomalous eastward current from the western Pacific converges with a westward current from the central-eastern Pacific near the dateline. This convergence drives the downwelling of warm surface water, creating a subsurface warm anomaly in the western-central Pacific (Fig. 6g). Collectively, these processes establish a characteristic vertical thermal structure marked by surface cooling overlying subsurface warming (Fig. 6g)—a signature of the slow oceanic response to persistent rainfall deficits. These dynamical responses are supported by prior idealized experiments25,48 and are validated in our GC-OGCM simulations (Supplementary Fig. 7).

Thus, under sustained 2-year rainfall deficits, the slow adjustment of oceanic processes causes prolonged equatorial cooling. Simultaneously, the rainfall-induced MLS anomalies facilitate the eastward propagation of successive cold signals, which collectively enhance the cooling in the second-year La Niña.

Methods
Observational and reanalysis datasets
To examine the MLS and rainfall response to multiyear La Niña events, a monthly mean comprehensive dataset derived from various methodologies is employed, including objective analyses and reanalysis data. The objective analysis data are from the Institute of Atmospheric Physics (IAP) global ocean temperature gridded product at a 1°×1° horizontal resolution (version 4)60, as well as from the EN.4.2.2 quality-controlled ocean dataset61. The reanalysis data were from the European Centre for Medium-Range Weather Forecasts Ocean Reanalysis System 5 (ORAS5)62, which provides three-dimensional potential temperature, salinity, velocity, and sea surface flux data. The sea surface temperature (SST) data were from the Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5)63. Owing to limitations in the observational record length of rainfall data from the Global Precipitation Climatology Project (GPCP)64, the study period spans from 1979 to 2022. Additionally, daily mean fields from the National Centers for Environmental Prediction-National Center for Atmospheric Research (NCEP–NCAR) reanalysis65 are used as the atmospheric fields to force the ocean model.

This study focuses on six multiyear La Niña events, specifically those occurring from 1983 to 1984, 1995 to 1996, 1998 to 2000, 2007 to 2008, 2010 to 2011, and 2020 to 2022; only the first 2 years of these events were included to maintain temporal consistency in the composite. Notably, the La Niña events of 2016 and 2017 are excluded from this analysis because of their relatively weaker amplitudes, as indicated by the Niño3.4 index (Fig. 1a). The long-term trends in all the data have been removed. In this study, the year when La Niña first develops is denoted year (0), the year when La Niña re-intensifies is denoted year (1), and the year when the multiyear La Niña decays is denoted year (2).

Ocean mixed-layer salinity budget
In this study, we employed an MLS budget derived from monthly mean ORAS5 reanalysis to assess the relative contributions of advection, vertical mixing, and freshwater flux (E-P) to the variations in salinity within the mixed layer. For the MLS budget,
∂S/∂t = −(E−P)S/h − u∂S/∂x − v∂S/∂y − w∂S/∂z + R (1)
where u, v, and S represent the mixed layer-averaged zonal current, meridional current, and salinity, respectively. w denotes the vertical velocity at the base of the mixed layer; E-P refers to the freshwater flux, E and P indicate evaporation and rainfall, respectively. R is the residual term, which accounts for contributions from unresolved processes such as horizontal and vertical diffusion; these processes cannot be directly estimated from the ocean reanalysis output. h indicates the monthly mean mixed layer depth (MLD), which is calculated with a surface-density-difference criterion of Δρ=0.03 kg m−366. A similar mixed-layer heat budget analysis has been performed in previous studies43. Furthermore, Eq.1 can be divided to emphasize the contributions of different ocean dynamics to the evolution of MLS as follows:
∂S0/∂t = −(E−P)S0/h − u0∂S̅/∂x − u′∂S0/∂x − u0∂S0/∂x − v0∂S̅/∂y − v′∂S0/∂y − v0∂S0/∂y − w0∂S̅/∂z − w′∂S0/∂z − w0∂S0/∂z + R (2)
where the overbars and primes denote the climatological mean and anomalous field, respectively.

Decomposing the individual contributions of temperature and salinity to density and stratification
We adopted a decomposition method to analyze the contributions of potential temperature and salinity to density and stratification37,67–69. The squared buoyancy frequency (N2) is calculated to represent stratification: N2= −g/ρ ∂ρ/∂z, where ρ and g indicate seawater density and gravitational acceleration, respectively. The individual contributions from potential temperature or salinity to N2 are calculated as follows:
N2_T = N2(T_inter, S_clim) − N2(T_clim, S_clim) (3)
N2_S = N2(T_inter, S_inter) − N2(T_inter, S_clim) ... (4)
where T_inter and S_inter indicate interannually varying potential temperature and salinity fields, respectively; T_clim and S_clim represent climatological potential temperature and salinity fields during 1980 and 2022, respectively. N2_T represents the contribution of interannual temperature anomalies to stratification variability, and N2_S analogously represents the contribution of salinity anomalies.

Ocean general circulation model and experiments
In this study, we used two ocean general circulation models (OGCMs) to investigate the impacts of rainfall and salinity on multiyear La Niña events. The models include a Gent-Cane (GC) OGCM45, which is based on sigma coordinates; and a Regional Ocean Modeling System (ROMS) model, which utilizes topography-following coordinates46. Both models successfully capture the spatiotemporal evolution characteristics of multiyear La Niña events. (Supplementary Fig. 8). Comparing these two models enhances the robustness of our analysis results. Detailed information about the specific configurations of the models and experiments are provided below (Supplementary Table 3).

GC-OGCM
The GC-OGCM is a reduced-gravity, primitive equation, sigma-coordinate model that encompasses the entire tropical Pacific (30°S-30°N, 120°E-76°W). This OGCM contains a mixed layer at the surface and 19 vertical layers beneath it according to sigma coordinates. A key characteristic of the model is that the mixed layer depth is prognostically calculated by a bulk mixed layer model developed by ref. 70. Additionally, the zonal resolution of the ocean model varies, with a resolution of 1° in the central basin, refining to 0.4° near the western and eastern boundaries. The meridional resolution is set to 1/3° near the equator and gradually increases to 2° poleward of 20° latitude. A sponge layer is implemented at the northern and southern boundaries from 20° to 30°. In this sponge layer, the temperature and salinity are gradually relaxed back to observed climatological fields. In the OGCM, the heat flux is explicitly calculated via an atmospheric mixed layer model, where only the wind speed, cloud cover, and solar radiation are prescribed71,72, consequently, latent and sensible heat fluxes are interactively determined by the surface ocean conditions in the OGCM. In the GC-OGCM, freshwater flux entering the ocean influences oceanic processes through two primary pathways: first, as a virtual salt flux; and second, by modifying the vertical velocity at the top layer via the natural boundary condition, which subsequently modulates the exchange of heat and momentum26,73. The model is initialized by the World Ocean Atlas (WOA) fields and then spin-up for 30 years forced by the NCEP/NCAR climatological fields (i.e., wind stress, radiation, and freshwater flux).

ROMS
ROMS is a free-surface, hydrostatic ocean model that solves primitive equations employs terrain-following (s-) coordinates46. The model domain encompasses the entire tropical Pacific and North Pacific regions (40°S–65°N, 110°E–70°W), with a horizontal resolution of 1/3°×1/3° to adequately resolve equatorial currents in the equatorial region, and it features 50-s-coordinate levels in the vertical direction69. Initial conditions are derived from the World Ocean Atlas 2009 (WOA09) fields and subsequently undergoes a 30-year spin-up period forced by the NCEP/NCAR climatological fields (i.e., wind stress, net heat flux, and freshwater flux)69. Following the spin-up, the model is forced by daily NCEP/NCAR atmospheric fields, from which surface fluxes are calculated online using bulk formulae.

Experiment I: Isolating the contributions of rainfall to MLS
In the control simulation (rain_inter), rainfall is prescribed as an interannually varying field. Both the GC-OGCM and ROMS are forced from 1979 to 2022 with daily NCEP/NCAR atmospheric data, from which surface fluxes are calculated online using bulk formulae. In the sensitivity experiment (rain_clim), rainfall is prescribed as its climatological mean for the period 1979–2022. The differences between the rain_inter and rain_clim experiments represent the oceanic response to interannually varying rainfall.

Given that the contribution of heat flux to MLS variability is relatively small over short timescales and acts primarily through its modulation of buoyancy flux and density (which subsequently influence mixing and circulation), the direct driver of MLS variability in the rain_clim experiment is ocean dynamical processes induced by wind forcing. Therefore, we did not explicitly distinguish the contributions of heat and momentum fluxes in this experimental framework.

Experiment II: Impacts of interannual rainfall anomalies on multiyear La Niña
To examine the impacts of rainfall-induced salinity changes on the ocean physical state, the GC-OGCM was forced with composite surface winds from NCEP/NCAR during multiyear La Niña events over a two-year period. A series of experiments were conducted as follows (Table S3):
In the climatological rainfall experiments, the observed monthly climatological rainfall (1980–2022) was prescribed to influence upper-ocean vertical velocity and salinity over the 2-year integration. In the first set of sensitivity experiments, the composite interannual rainfall anomalies during multiyear La Niña events were doubled and superimposed onto the climatological rainfall field (hereafter referred to as the Rainfall X2 experiment). The difference between the Rainfall X2 and climatological rainfall experiments characterizes the effect of amplified rainfall anomalies on ocean salinity and temperature. The doubling of rainfall anomalies was necessary because the GC-OGCM, due to its specific configuration, simulates weaker interannual MLS variability than observed. In this setup, rainfall influences both mixed-layer salinity and velocity (via the natural boundary condition)26,74; thus, the experiment confirms the effect of rainfall on both MLS and SST. A parallel experiment using the original, undoubled composite rainfall anomalies (Rainfall X1) was also conducted to test sensitivity.

In the second set of sensitivity experiments, a constant rainfall anomaly of –9 mm day⁻¹ was artificially imposed on the climatological rainfall field in the western-central equatorial Pacific (140°E–180°, 5°S–5°N) from June(0) to December(1) of the multiyear La Niña event (referred to as the constant rainfall experiment). In this configuration, the rainfall anomaly was designed to influence ocean salinity exclusively, without modifying the initial vertical velocity (i.e., the natural boundary condition). This was achieved by disabling the impact of rainfall anomalies on vertical velocity within the model code. The region and magnitude of the rainfall anomaly were chosen to closely reproduce the observed MLS anomalies. The results show that the simulated salinity response is consistent with observations, though some differences remain—as analyzed in the Results section.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Multi-year La Niña events occur frequently, and have a significant impact on global weather and climate. The present study examined multi-year La Niña events, which are self-sustained by changes in rainfall and salinity. The authors revealed the relationships between rainfall, salinity and the SST in both the central and eastern Pacific based on extensive data analysis and model experiments. I particularly appreciated the authors’ efforts on model simulations. They conducted two sets of ocean general circulation model experiments to illustrate the delayed effect of rainfall on salinity, and a series of well-designed sensitivity experiments to examine the effect of rainfall on SST via salinity during multi-year La Niña events. The model results provide valuable evidence to support effects of rainfall and salinity on multi-year La Niña events. The paper is well-written, and the methods are convincing. Although the role of rainfall-salinity feedback in multi-year La Niña events is high novel, some revisions are needed to clarify some issues before its publication.\nSpecific comments:\n1. The abstract requires a clearer description of how salinity responds to multi-year La Niña events. It should also specify the relative contributions of rainfall and oceanic dynamic processes during multi-year La Niña events. The explanation of the slow response process to rainfall is overly simplistic and requires a more comprehensive introduction. Furthermore, the implications presented in the abstract need to be strengthened to convey the impact of this study effectively.\n2. Model agreement and discrepancies in MLS simulation.\nI note that the simulated MLS variability in the GC model is weaker than that in ORAS5 reanalysis or IAP observational data. While the authors have offered plausible explanations for this discrepancy, the potential implications of this model bias for the subsequent sensitivity experiments require further clarification. For instance, how might this amplitude difference affect the intensity of rainfall-salinity feedback? Regarding the ROMS simulation, although the amplitude of the modeled MLS matches the observations well, notable differences exist in the spatial patterns. Specific reasons for these spatial discrepancies should be addressed.\n3. In \"Asynchronous effects of… on MLS\" section, some of the texts overlaps with previous section and should be condensed. While the significant role of zonal advection is acknowledged, a more detailed discussion of the impact of vertical mixing processes is needed. A more detailed analysis of the balance between rainfall contribution and horizontal advection in Figure 2a would be beneficial. Specifically, decomposing horizontal advection into its zonal and meridional components during the NDJF period could reveal important underlying physical processes related to oceanic adjustments. The current figure layout, while informative, could be enhanced through reorganization to better highlight these causal relationships.\n4. The paper used ORAS5 data for the diagnostic analysis. Why not consider using IAP4 or EN4 based on the observations? What differences might arise in the results? Additionally, I am wondering what role salinity plays in stratification during a La Niña event in a single year, and how this compares with events spanning multiple years (double dip or triple dip)?\n5. A more comprehensive discussion of the impact of rainfall and ocean salinity on ENSO diversity will benefit the manuscript. In particular, the influence of rainfall on physical processes after it enters the ocean requires further clarification, despite being mentioned briefly in the Methods section.\n6. I am unsure about the results presented in Extended Data Figure 3. The discussion is rather simplistic and lacks depth. To improve clarity, I recommend incorporating this figure into the main text.\n7. Lines 76-86: The paper emphasizes that rainfall anomalies caused by multi-year La Niña events can have a lasting impact on the ocean, which is further intensified due to the changes in ocean salinity. However, the paragraph does not clearly explain their cumulative effect, which I believe is essential for readers to understand.\n8. Fig.S1: GODAS shows relatively weak interannual anomalies for MLS. What might be causing of this? Please investigate this further and confirm the cause.\n9. Line 156: I suggest rewriting this as “Delayed response of MLS to rainfall”.\n10. Figure 7: Due to the complexity of how rainfall enters the ocean, a schematic diagram is an effective way to illustrate the relevant processes. The current diagram is overly complicated, and not mentioned in main texts. It would be better to focus on the key processes, while omitting relevant details.", "reviewer_2": "Comments on “Rainfall sustains multiyear La Niña\" (NCOMMS-25-41284) by Tian et al.\nThis paper explored the effects of rainfall in sustaining multi-year La Nina by conducting OGCM experiments. It was found that the positive salinity anomalies in the western-central Pacific can enhance La Nina amplitude significantly. Considering the multi-year La Nina becomes more frequent, physics behind need to be explored. This study made an important contribution to the topic. I suggest it be accepted after some revisions/clarifications.\nMajor point:\nThe authors are suggested to include some discussions about robustness of your finding. For example, how much are your quantitative estimate of 14% and 32% amplitude enhancement by the rainfall reduction depending on the model you are using? As seen from Extended Data Figure 2 and Table 1, GC-OGCM presents a significant deficit in simulating multi-year La Nina, particularly in simulating their first-year peak. How does ROMS perform in simulating multi-year La Nina. I am actually curious about how multi-year La Nina events were simulated in the control simulation raininter (L477; Sect c.)\nOther points:\nL113/Fig S1/L394-395: As a very old system, GODAS does not assimilate real salinity data. It makes no sense to include it in such diagnostics.\nL123-124: While previous budget analyses did suggest the critical roles of the two terms, your Fig. 2a seems not. Can you give some explanations?\nL167: Other than saying “systematic biases”, it is more informative to specify what biases you are referring to.\nL200: Figs 4a and 4b=>4b and 4c;\nL204: winter=>summer; Fig. 4c=>4a\nL209: after “instability”, add “in the western Pacific”\nL247: Can you explain how you derive the numbers? In Extended Data Table 1, it is useful to include the reference information in parallel to the observed values.\nL414: Are your budget analyses based on monthly or daily outputs?\nPlease explain how the freshwater flux (evaporation minus precipitation) is represented in your OGCMs? By natural boundary condition (Huang 1993)?\nL484/Sect. d: Ensemble information is needed for the experiments.\nFig1: Which datasets are used for the analyses?"}, "author_response": "Response to Reviewers\nReviewer #1:\nMulti-year La Niña events occur frequently, and have a significant impact on global\nweather and climate. The present study examined multi-year La Niña events, which are\nself-sustained by changes in rainfall and salinity. The authors revealed the relationships\nbetween rainfall, salinity and the SST in both the central and eastern Pacific based on\nextensive data analysis and model experiments. I particularly appreciated the authors’\nefforts on model simulations. They conducted two sets of ocean general circulation\nmodel experiments to illustrate the delayed effect of rainfall on salinity, and a series of\nwell-designed sensitivity experiments to examine the effect of rainfall on SST via\nsalinity during multi-year La Niña events. The model results provide valuable evidence\nto support effects of rainfall and salinity on multi-year La Niña events. The paper is\nwell-written, and the methods are convincing. Although the role of rainfall-salinity\nfeedback in multi-year La Niña events is high novel, some revisions are needed to\nclarify some issues before its publication.\n=> Thank you very much for your comments and suggestions. Please find our point-\nby-point responses (in blue font) below, as well as the revisions (in red font) in the\nMain Text, all addressing your feedback. In addition, we have thoroughly revised and\npolished the entire text for clarity and consistency.\nSpecific comments:\n1. The abstract requires a clearer description of how salinity responds to multi-year La\ndynamic processes during multi-year La Niña events. The explanation of the slow\nresponse process to rainfall is overly simplistic and requires a more comprehensive\nintroduction. Furthermore, the implications presented in the abstract need to be\nstrengthened to convey the impact of this study effectively.\n Thank you very much for your suggestions. We have revised the abstract: we added\nan introduction to the specific processes affecting mixed-layer salinity (MLS),\nenriched the description of the processes in the second year of multi-year La Niña\nevents, and expanded the discussion on implications. Abstract has been modified as\nfollows, where the words in red denote the revisions in Lines 17-35.\n“Multiyear La Niña events are becoming more frequent, as shown by historical\nobservational records and climate models. These prolonged events exert cumulative\nimpacts on tropical Pacific rainfall variability, which subsequently affects mixed-layer\nsalinity (MLS) in the region. However, whether and how cumulative rainfall anomalies\ninfluence multiyear La Niña events through a salinity effect remains poorly understood.\nHere, using observations and ocean model experiments, we show that the positive MLS\nanomalies in the western–central equatorial Pacific—driven initially by oceanic\ndynamical processes and later by rainfall deficits during multiyear La Niña—act to\nreinforce these prolonged cooling events. Quantitatively, rainfall reduction enhances La\nNiña amplitude by 14% in the first year and 32% in the second year. Initially, MLS\nanomalies in the western–central equatorial Pacific trigger equatorial Kelvin wave\n2\nadjustments, inducing rapid surface cooling in the eastern equatorial Pacific.\nSubsequently, slower ocean circulation responses lead to basin-wide equatorial Pacific\ncooling. The rapid oceanic response to rainfall anomalies, superimposed on the slow\ncooling response, generates a cumulative positive feedback from rainfall on the second-\nyear La Niña. These findings reveal rainfall–salinity feedbacks as a key mechanism\nsustaining multiyear La Niña events, with broad implications for ENSO prediction and\nthe global hydrological cycle.”\n2. Model agreement and discrepancies in MLS simulation.\nI note that the simulated MLS variability in the GC model is weaker than that in ORAS5\nreanalysis or IAP observational data. While the authors have offered plausible\nexplanations for this discrepancy, the potential implications of this model bias for the\nsubsequent sensitivity experiments require further clarification. For instance, how\nmight this amplitude difference affect the intensity of rainfall-salinity feedback?\nRegarding the ROMS simulation, although the amplitude of the modeled MLS matches\nthe observations well, notable differences exist in the spatial patterns. Specific reasons\nfor these spatial discrepancies should be addressed.\n 1) We sincerely appreciate your suggestions. In this revision, we have added a series\nof sensitivity experiments (24 experiments) to verify the impacts of spatial and\ntemporal differences in the rainfall pattern imposed in the GC-OGCM on multiyear\nLa Niña events. Overall, changes in the rainfall pattern do not alter the sustaining\neffect of the rainfall-salinity feedback on multiyear La Niña. Specifically, the\nexperiments are divided into three types: applying rainfall anomalies in the first-\nyear La Niña (1yr exp), applying rainfall anomalies in the second-year La Niña\n(2yr exp), and applying rainfall anomalies throughout the entire period of multiyear\nLa Niña (all exp). The rainfall anomaly is set at a fixed value of -9 mm/day to ensure\ncomparability with the observed MLS amplitude. On this basis, we have performed\na zonal shift of the region where rainfall is applied to test the regional dependence\nof the rainfall-salinity feedback. Totally, eight sets of experiments are designed,\neach centered in the range of 150°E–140°W, with a zonal extent of 40°. For instance,\nthe experiment centered at 150°E covers a spatial range of 130°E–170°E.\n2) Firstly, the experimental results show that in the 1yr exp, the impact of rainfall\non MLS in the Niño4 region is more significant in the first-year La Niña (Figure\nR1); however, a distinctly positive anomaly persists in the second year, indicating\nthat the initial rainfall can still affect MLS in the western-central equatorial Pacific\nin the subsequent year. Additionally, in the 2yr exp, the impact of rainfall on MLS\nin the second year is stronger than that of the 1yr exp on MLS in the first year,\nsuggesting that the background field in the second year of multi-year La Niña is\nmore conducive to amplifying salinity changes induced by rainfall. When\nconsidering the all exp, the impact of rainfall on MLS in the second year represents\nthe superposition of the effects from the aforementioned two experiments, which\nfurther confirms the existence of a cumulative effect of rainfall on MLS during\nmultiyear La Niña. Secondly, from the perspective of spatial differences, the mixed\nlayer in the eastern equatorial Pacific is relatively shallow, making it prone to a\nstronger MLS response to rainfall anomalies; moreover, under the influence of\n4\nwestward currents such as the South Equatorial Current (SEC), applying rainfall in\nthe central-eastern Pacific can induce more significant salinity anomalies in the\nNiño4 region. As you can see, the MLS anomalies simulated by the GC-OGCM are\nmainly concentrated in the observed western-central equatorial Pacific; thus, the\nbiases introduced by the simulation are relatively small.\nFigure R1 | Responses of local MLS and SST to varying rainfall patterns. In a,\na -9 mm/day rainfall anomaly is imposed in the GC-OGCM to influence MLS, with\n5\nthe anomaly center location varied (e.g., a center at 150°E corresponds to a zonal\nregion of 130°E–170°E [40° in width] and a meridional range of 5°S–5°N). “1yr\nexp” denotes experiments where the rainfall anomaly is imposed during the first-\nyear La Niña, spanning June (0) to February (1); “2yr exp” denotes experiments\nwhere the anomaly is imposed during the second-year La Niña, spanning March (1)\nto December (1); “all exp” denotes experiments where the anomaly is imposed\nthroughout the entire multiyear La Niña period, spanning June (0) to December (1).\nTwo peak phases—December (0)–January–February (1) [D(0)JF(1)] and December\n(1)–January–February (2) [D(1)JF(2)]—are selected for analysis. In b, Niño 3.4\nindexes from the different experiments are presented.\n 3) Regarding the ROMS simulations, we also compared the differences in sea\nsurface temperature (SST) simulations during multiyear La Niña events (Figure R2).\nThe results show that both the GC-OGCM and ROMS can well reproduce the\nspatial anomaly distribution of SST during multi-year La Niña; however, both\nmodels underestimate the intensity of SST cold anomalies in the eastern equatorial\nPacific. This, on the one hand, is a common issue in current ocean models, mainly\nattributed to the weak simulated vertical mixing processes and zonal currents. Such\nweaknesses may lead to reduced zonal advection of salinity and temperature in the\neastern equatorial Pacific, ultimately resulting in underestimated SST simulations\non the one hand and overestimated MLS simulations on the other.\n2°S-2°N for the multi-year La Niña events, derived from (a) observation from OBS\n(ERSSTv5), and (b) OGCM simulation, (c) ROMS during 1980-2022. Stippling\nrepresents values above the 90% confidence level by a two-tailed t-test.\n 4) Some texts are added as follows (red fonts in Lines 163-166; 259-265 from\nMain Text):\n“However, due to the weaker vertical mixing and entrainment associated with wind\nstirring in the GC-OGCM47, the simulated amplitude of interannual MLS anomalies\nis subsequently underestimated.”\n“Furthermore, additional experiments with rainfall anomalies applied at different\nlongitudes confirm that the amplitudes of the SST and MLS responses during\nmultiyear La Niña events are generally consistent across different spatial\ndistributions of rainfall (Fig. S6). These experiments collectively demonstrate that\nthe effects of negative rainfall anomalies from the first and subsequent years\nsuperimpose, contributing to the more pronounced cooling in the second-year La\nNiña.”\n3. In \"Asynchronous effects of… on MLS\" section, some of the texts overlaps with\nprevious section and should be condensed. While the significant role of zonal advection\nis acknowledged, a more detailed discussion of the impact of vertical mixing processes\nis needed. A more detailed analysis of the balance between rainfall contribution and\nhorizontal advection in Figure 2a would be beneficial. Specifically, decomposing\nhorizontal advection into its zonal and meridional components during the NDJF period\ncould reveal important underlying physical processes related to oceanic adjustments.\nThe current figure layout, while informative, could be enhanced through reorganization\nto better highlight these causal relationships.\n=> Thank you for your suggestions. We have significantly condensed and integrated\nthe first two sections to highlight the significant contributions of rainfall and ocean\ndynamical processes during different phases of multiyear La Niña. Regarding vertical\nmixing processes, we consider this a critical term influencing the evolution of MLS,\nwhich is particularly crucial in the early stage of the first-year La Niña. We believe this\nis somewhat related to the deepening of the mixed layer during the initial phase of La\nNiña (Figure R3) due to enhanced vertical mixing and entrainment. However, as you\nhave observed, we performed an offline budget analysis using ORAS5 data, so the\nvertical mixing term can only be treated as a residual. Thus, the role of this term may\nnot be overinterpreted. A more detailed investigation into the combined effects of\n8\nvertical mixing and zonal advection requires further clarification in future studies.\nConcerning meridional advection, we also performed an analysis and found that its\nposition in the ITCZ does contribute somewhat to MLS, but it remains a minor term\nthroughout the entire multi-year La Niña period (Figure R4). Meanwhile, we have\nfollowed your suggestion and redrawn Fig. 2 to better present the logic and timeline.\nFigure R3| Contribution from vertical mixing and other processes (Units: psu month-1)\nas the residual terms to MLS tendency in MAMJ (0), ND(0)JF(1), and ASON (1),\nrespectively.\nFigure R4| Meridional advection of the mean salinity by the anomalous current ( )\n′\n(Units: psu month-1) in MAMJ (0), ND(0)JF(1), and ASON (1), respectively. 𝑣 𝑆̅\n4. The paper used ORAS5 data for the diagnostic analysis. Why not consider using\nIAP4 or EN4 based on the observations? What differences might arise in the results?\nevent in a single year, and how this compares with events spanning multiple years\n(double dip or triple dip)?\n=> Thank you for your suggestions. We have redrawn Fig. 4 using IAP data. It is found\nthat the IAP data exhibit an almost identical spatiotemporal evolution to that of ORAS5;\nthe only difference is that the impact of salinity on N² at the base of the mixed layer is\nweaker in IAP data than in ORAS5. Nevertheless, the IAP data can basically reproduce\nthe results of ORAS5 and also illustrate the key role of salinity in determining the\nnegative anomalies of N² at the base of the mixed layer.\nRegarding the comparison with single-year La Niña events, we consider this an\nexcellent suggestion. However, since most La Niña events since 1980 have occurred\nconsecutively over two years, there are only two isolated single-year La Niña events\n(i.e., 1988-1989 and 2005-2006). Therefore, to better demonstrate the impact of salinity\non N² during the first year of La Niña, we conducted a composite analysis of the first\nyear of all La Niña events. The results show that the potential temperature begins to\ndecrease from the winter of the La Niña year, and the cold temperature anomalies\ndisappear by the summer of the following year. Surprisingly, the positive salinity\nanomalies can persist until the winter of the second year, when the SST has already\nreturned to normal conditions. These prolonged positive salinity anomalies further\naffect the stratification, which can promote vertical mixing and thereby contribute to\nthe occurrence of the second-year La Niña. This indirectly confirms the sustaining\n10\neffect of salinity on multi-year La Niña events.\nFigure R5| Diagnostic analysis for salinity effect on the upper ocean. a, Time‒depth\n(0–200 m) section for density anomalies averaged over the western equatorial Pacific\nduring multiyear La Niña. b-d, Same as a but for temperature, salinity, and squared\nbuoyancy frequency (N2), respectively. e-f, N2 changes due to interannual temperature\nanomalies and salinity anomalies, respectively. The green (black) contours in a–f\nrepresents the MLD from the climatological mean and composite multiyear La Niña.\nThe potential temperature and salinity data are from the IAP datasets. Stippling\nrepresents values above the 90% confidence level according to a two-tailed t test.\nfirst-year La Nina of all La Nina events. a, Time‒depth (0–200 m) section for density\nanomalies averaged over the western equatorial Pacific during first-year La Niña. b-d,\nSame as a but for temperature, salinity, and squared buoyancy frequency (N2),\nrespectively. e-f, N2 changes due to interannual temperature anomalies and salinity\nanomalies, respectively. The green (black) contours in a–f represents the MLD from the\nclimatological mean and composite first-year La Niña. The potential temperature and\nsalinity data are from the ORAS5 datasets. Stippling represents values above the 90%\nconfidence level according to a two-tailed t test.\n5. A more comprehensive discussion of the impact of rainfall and ocean salinity on\nENSO diversity will benefit the manuscript. In particular, the influence of rainfall on\nphysical processes after it enters the ocean requires further clarification, despite being\nmentioned briefly in the Methods section.\n Thank you very much for your suggestions. As previous study shown, an increase\n12\nin MLS leads to a rise in upper-layer density, which in turn further induces\nstratification instability and enhances vertical mixing. This process thereby\nfacilitates the entrainment and mixing of subsurface cold water into the surface layer,\nultimately resulting in cooling of the surface layer. (1) There have been several\nstudies on the impact of salinity changes induced by rainfall or freshwater flux on\nENSO diversity, which mainly based on using the GC-OGCM as adopted in this\nstudy. Therefore, we have also included some introductions to this issue in the\nIntroduction section. (2) In Regarding how rainfall affects oceanic processes after\nentering the ocean, we have also provided some descriptions, as detailed below\n(Lines in 56-59; 429-432).\n“Nevertheless, during multiyear La Niña events, the equatorial ocean undergoes\nsustained cooling for over two years—a timescale that implies slow ocean\ndynamical processes (e.g., vertical mixing, wind-driven large-scale circulation, and\nrainfall-induced salinity changes) may play a pivotal role.\n“In the GC-OGCM, freshwater flux entering the ocean influences oceanic processes\nthrough two primary pathways: first, as a virtual salt flux; and second, by modifying\nthe vertical velocity at the top layer via the natural boundary condition, which\nsubsequently modulates the exchange of heat and momentum26,75.”\n6. I am unsure about the results presented in Extended Data Figure 3. The discussion\nis rather simplistic and lacks depth. To improve clarity, I recommend incorporating\nthis figure into the main text.\n13\n We sincerely appreciate your suggestions. We indeed believe that Extended Figure\n3 is helpful for understanding the cumulative effect of rainfall on multiyear La Niña\nevents. Meanwhile, similar to our response to your second question, we have\nconducted additional sets of experiments to investigate the responses of multi-year\nLa Niña to the spatial and temporal distributions of rainfall. Therefore, we have\nincorporated some additional discussions into the main text. However, as more\nexperimental results might make the paper excessively lengthy and less focused on\nkey physical processes, we have moved some figures to the supplementary\nmaterials. We hope you understand our intentions. Some texts are added as follows\n(Lines 259-265; 315-318 in the Main Text):\n“Furthermore, additional experiments with rainfall anomalies applied at different\nlongitudes confirm that the amplitudes of the SST and MLS responses during\nmultiyear La Niña events are generally consistent across different spatial\ndistributions of rainfall (Fig. S6). These experiments collectively demonstrate that\nthe effects of negative rainfall anomalies from the first and subsequent years\nsuperimpose, contributing to the more pronounced cooling in the second-year La\nNiña.\n…Thus, under sustained two-year rainfall deficits, the slow adjustment of oceanic\nprocesses causes prolonged equatorial cooling. Simultaneously, the rainfall-induced\nMLS anomalies facilitate the eastward propagation of successive cold signals,\nwhich collectively enhance the cooling in the second-year La Niña.”\n14\n7. Lines 76-86: The paper emphasizes that rainfall anomalies caused by multi-year\nLa Niña events can have a lasting impact on the ocean, which is further intensified\ndue to the changes in ocean salinity. However, the paragraph does not clearly\nexplain their cumulative effect, which I believe is essential for readers to understand.\n Thank you for your suggestions. We have revised this section to more accurately\nconvey the cumulative effect (Lines in 85-97).\n“Therefore, it is hypothesized that since rainfall-induced MLS anomalies can persist\nfor over two years during multi-year La Niña events, they may exert a sustained\ninfluence on the upper ocean, which, in turn, could sustain the cooling. The\nunderlying processes are plausible: the absence of freshwater input from negative\nanomalies in rainfall reduces upper-ocean stratification, thereby enhancing vertical\nmixing and entrainment, deepening the mixed layer, and modulating large-scale\nocean circulation and equatorial waves. These processes collectively regulate SST\nevolution over extended timescales in a manner that favors the persistence of La\nNiña. However, the extent to which these rainfall-driven MLS variations govern the\nprogression of multi-year La Niña events remains undetermined. In this study, we\naddress this issue using observational data and numerical simulations,\ndemonstrating that multi-year La Niña events exhibit a self-amplifying nature\nthrough a cumulative rainfall–salinity feedback mechanism.”\n8. Fig.S1: GODAS shows relatively weak interannual anomalies for MLS.\nWhat\nmight be causing of this? Please investigate this further and confirm the cause.\n15\n Thank you very much for your suggestions. Since GODAS does not assimilate\nobserved salinity data, the interannual variability of MLS simulated by GODAS is\nabnormally weak. Therefore, we have removed the content related to GODAS in\nthis revision. Thank you.\n9. Line 156: I suggest rewriting this as “Delayed response of MLS to rainfall”.\n=> Corrected as suggested. Thank you\n10. Figure 7: Due to the complexity of how rainfall enters the ocean, a schematic\ndiagram is an effective way to illustrate the relevant processes. The current diagram is\noverly complicated, and not mentioned in main texts. It would be better to focus on the\nkey processes, while omitting relevant details.\n=> Thank you for your suggestions. We have redrawn the schematic diagram.\nFigure 7 | Schematic diagram of rainfall–salinity feedbacks sustaining multiyear\nLa Niña. a, Initial phase: In the western-central equatorial Pacific, positive mixed-layer\nsalinity (MLS) anomalies are initially driven by zonal advection and vertical mixing.\nThese anomalies deepen the mixed layer, exciting a baroclinic response with alternating\n16\nsubsurface warm/cold layers. An associated cold anomaly propagates eastward via\nequatorial wave dynamics, inducing remote surface cooling in the eastern basin. b,\nMature phase: Persistent rainfall reduction now dominates the maintenance and\namplification of MLS anomalies. These sustained salinity anomalies drive large-scale,\nslow oceanic adjustments—including modified circulation and enhanced vertical\nmixing—that reinforce basin-wide cooling and prolong the La Niña state.\nComments on “Rainfall sustains multiyear La Niña\" (NCOMMS-25-41284) by Tian et\nal.\nThis paper explored the effects of rainfall in sustaining multi-year La Nina by\nconducting OGCM experiments. It was found that the positive salinity anomalies in the\nwestern-central Pacific can enhance La Nina amplitude significantly. Considering the\nmulti-year La Nina becomes more frequent, physics behind need to be explored. This\nstudy made an important contribution to the topic. I suggest it be accepted after some\nrevisions/clarifications.\nThank you very much for your comments and suggestions. Please find our point-by-\npoint responses (in blue font) below addressing your feedback. In addition, we have\nthoroughly revised and polished the entire text. Please refer to the Main Text for\ndetailed revisions (in red font).\nMajor point:\nThe authors are suggested to include some discussions about robustness of your finding.\nFor example, how much are your quantitative estimate of 14% and 32% amplitude\nenhancement by the rainfall reduction depending on the model you are using? As seen\nfrom Extended Data Figure 2 and Table 1, GC-OGCM presents a significant deficit in\nsimulating multi-year La Nina, particularly in simulating their first-year peak. How\ndoes ROMS perform in simulating multi-year La Nina. I am actually curious about how\nmulti-year La Nina events were simulated in the control simulation raininter (L477;\nSect c.)\n18\n We greatly appreciate your suggestions. Regarding the issue of result robustness as\nyou mentioned, we will respond to each point one by one. First, concerning the\nmodel dependence in the contribution of rainfall to the amplitude of multi-year La\nNiña events in the OGCM results, we have conducted in-depth discussions. We\nperformed a series of sensitivity experiments to verify the impacts of imposing\nrainfall at different periods and locations on multi-year La Niña events (Figure R1).\nSecond, we summarize existing studies on the impacts of incorporating rainfall\nanomalies into different ocean models or coupled models on the mean state of the\ntropical Pacific and ENSO (Table R1). Third, we conducted experiments based on\nthe Regional Ocean Modeling System (ROMS). The results indicate that the\nmagnitudes of the rainfall-salinity effects in the GC-OGCM and ROMS (~0.1-0.2°\nC) are basically comparable to that in other models (Table R1 and Figure R2).\nFinally, we evaluated the ability of ROMS and GC-OGCM to simulate multi-year\nLa Niña events, and found that both models can basically reproduce the observed\nspatial patterns (Figure R3). Overall, we have added Table S1 and Figures S3–S4\nto evaluate the model's simulation capability and the sensitivity of the experimental\nresults. Corresponding descriptions have been incorporated into the main text. The\nspecific responses are as follows.\n1. Robust effect of rainfall on multiyear La Nina in GC-OGCM experiments\nThe experiments are divided into three types: applying rainfall anomalies in the\nfirst-year La Niña (1yr exp), applying rainfall anomalies in the second-year La Niña\n19\n(2yr exp), and applying rainfall anomalies throughout the entire period of multi-year\nLa Niña (all exp). The rainfall anomaly is set at a fixed value of -9 mm/day to ensure\ncomparability with the observed MLS amplitude, such that the impact of rainfall via\nthe salinity pathway is reasonably represented. Nevertheless, these idealized\nexperiments only emphasize the sensitivity of the rainfall effect to location and timing.\nOn this basis, we have performed a zonal shift of the region where rainfall is applied to\ntest the regional dependence of the rainfall feedback. Overall, eight sets of experiments\nare designed, each centered in the range of 150°E–140°W, with a zonal extent of 40°.\nFor instance, the experiment centered at 150°E covers a spatial range of 130°E–170°E.\nFirstly, the experimental results show that in the 1yr exp, the impact of rainfall on\nMLS in the Niño4 region is more significant in the first-year La Niña (Figure R1);\nhowever, a distinctly positive anomaly persists in the second year, indicating that the\ninitial rainfall can still affect MLS in the western-central equatorial Pacific in the\nsubsequent year. Additionally, in the 2yr exp, the impact of rainfall on MLS in the\nsecond year is stronger than that of the 1yr exp on MLS in the first year, suggesting that\nthe background field in the second year of multi-year La Niña is more conducive to\namplifying salinity changes induced by rainfall. When considering the all exp, the\nimpact of rainfall on MLS in the second year represents the superposition of the effects\nfrom the aforementioned two experiments, which further confirms the existence of a\ncumulative effect of rainfall on MLS during multiyear La Niña. Secondly, from the\nperspective of spatial differences, the mixed layer in the eastern equatorial Pacific is\nrelatively shallow, making it prone to a stronger MLS response to rainfall anomalies;\n20\nmoreover, under the influence of westward currents such as the South Equatorial\nCurrent (SEC), applying rainfall in the central-eastern Pacific can induce more\nsignificant salinity anomalies in the Niño4 region. As you can see, the MLS anomalies\nsimulated by the GC-OGCM are mainly concentrated in the observed western-central\nequatorial Pacific; thus, the biases introduced by the simulation are relatively small.\nFigure R1 | Responses of MLS and SST to varying rainfall patterns. In a, a -9\nmm/day rainfall anomaly is imposed in the GC-OGCM to influence MLS, with the\n21\nanomaly center location varied (e.g., a center at 150°E corresponds to a zonal region of\n130°E–170°E [40° in width] and a meridional range of 5°S–5°N). “1yr exp” denotes\nexperiments where the rainfall anomaly is imposed during the first-year La Niña,\nspanning June (0) to February (1); “2yr exp” denotes experiments where the anomaly\nis imposed during the second-year La Niña, spanning March (1) to December (1);\n“all exp” denotes experiments where the anomaly is imposed throughout the entire\nmultiyear La Niña period, spanning June (0) to December (1). Two peak phases—\nDecember (0)–January–February (1) [D(0)JF(1)] and December (1)–January–February\n(2) [D(1)JF(2)]—are selected for analysis. In b, Niño 3.4 indexes from the different\nexperiments are presented.\n2. Model dependency\nRegarding the potential model dependency issue arising from our use of a single\nmodel, since we are currently unable to perform unified experiments with multiple\nocean models, we have summarized as comprehensively as possible the results from\nprevious literature on the impacts of freshwater flux or rainfall on the tropical Pacific\nmean state and ENSO. Specifically, most studies have confirmed that when negative\nrainfall anomalies with an amplitude comparable to interannual variability are imposed\non the western tropical Pacific, they can lead to a cooling of more than -0.3°C in the\ncentral-western equatorial Pacific (Kang et al., 2017; Wu et al., 2010)—this range is\nconsistent with the SST changes observed in our current experiment. Alternatively,\nwhen interannually varying rainfall anomalies are incorporated into multiple ocean\n22\nmodels, it is found that they can increase the ENSO amplitude by 10-15%, a magnitude\nthat aligns with the observational results in this study (Guan et al., 2022, 2023; Maes et\nal., 2002; Vialard et al., 2002; Yang et al., 1999; Zhu et al., 2014).\nFurthermore, we have also conducted experiments similar to the GC-OGCM using\nthe ROMS. The results show that introducing negative rainfall anomalies in ROMS\nleads to a more pronounced cooling of the second-year La Niña (~0.1°C), though the\nintensity is lower than that in the GC-OGCM. The cooling effect of rainfall on SST in\nthe ROMS is relatively weak, which may be attributed to the fact that its atmospheric\nforcing field is fully prescribed. In contrast, in the GC-OGCM, the key components of\nnet heat flux (i.e., latent heat and sensible heat) are determined by an atmospheric mixed\nlayer model (Refs. 9-10); thus, partial atmospheric feedbacks have been represented in\nGC-OGCM. Another factor is that the uppermost layer in GC-OGCM is the mixed layer,\nwhich results in the cooling effect of rainfall being more readily reflected in SST.\nWhereas in ROMS, rainfall is only added to the topmost layer (~5 m), so the impact of\nrainfall needs to propagate downward gradually through vertical diffusion. Therefore,\ncombined with Table R1, we consider that our current conclusion—i.e., rainfall can\nmaintain approximately 14-32% of the SST cold anomaly amplitude during multi-year\nLa Niña events—is generally reasonable. However, in the discussion section, we also\npropose that more ocean models or climate models should be used to conduct the same\nnumerical experiments to further consolidate this conclusion.\nTable R1 |Summary of previous studies on numerical experiments regarding the impact\n23\nof rainfall on the tropical Pacific mean state and ENSO\nReference Model used Method for added Impacts or results\nrainfall\nWu et al., Fast Ocean Atmosphere Model -3.456 mm day-1 in 10°S–10°N, -0.3°C in the central-western tropical\n(2010) (FOAM); Ocean Model-MOM3 60°E–180° . Pacific\nVialard et A tropical Pacific version of the Interannual rainfall effect is 0.5°C warming in the equatorial ocean\nal., 2002 Laboratoire d’Oce´a- nographie investigated compared to the during 1997-1998 El Nino\nDYnamique et de Climatologie climatological rainfall effect.\n(LODYC) Ocean General\nCirculation Model\nMaes et al, Coupled model at the Laboratoire Modified N2 change with or In the absence of salinity stratification,\n2002 d’Oce´anographie Dynamique et without salinity effect in the slightly cooler SST (~0.5°) result in\nde Cli- matologie (LODYC) warm pool a reduced El Nino or a return to the\nmean seasonal cycle of the model.\nZhu et al., NCEP Climate Forecast System, The experimental predictions Tendency to decay of the subsurface\n2014 version 2 (CFSv2). are produced for up to 12 cold condition during the spring and\nmonths lead time from April to early summer 2007 was interrupted by\nthe following March, referred to positive salinity anomalies in the\nas the CTL experiment. In the upper central Pacific\nsensitivity, the initial salinity\nstates in experiment CTL were\nreplaced with the ORAS4\nsalinity climatology.\nKang et al., MOM5 an extended elliptical zone SST increases along the equator,\n2017 along the equator with the especially near the cold tongue and\nmaximum values (+8.0 around 170°W (approximately 0.2 °C)\nmm/day) located on the\nInternational Date Line\nYang et al., A quasi-isopycnal general Interannual rainfall effect is The magnitude of this gain in SST\n1999 circulation model as described in investigated. caused by the enhanced precipitation\nSchopf and Loughe (1995). is about 10-15% of the total SST\nincrease observed during the winter,\nGuan et al., GC-OGCM in this study Interannual rainfall effect with The distinct zonal structures of SSS\n2022;2023 mimicking observed salinity anomalies facilitate stronger EPEN\nvariability . than the CPEN, increasing their\ndifference in intensity by about 10%\ndue to the interannual anomalies in the rainfall in ROMS, respectively.\n3. Model Biases\nRegarding the issue of GC-OGCM simulation bias in Extended Figure 2 (Fig.S4\nin current version) that you mentioned, we sincerely apologize for failing to clearly\nindicate in the figure that each line graph in the time dimension represents the difference\nbetween the sensitivity experiment (e.g., Rainfall X2) and the control experiment (with\nrainfall set to the climatological mean state). This oversight led to a misunderstanding\nbetween the simulated SST changes and the observed amplitude. To address this, we\nhave now labeled on the axes which data series correspond to the differences between\nrespective experiments and which correspond to the actual SST amplitude during La\nNiña events.\nAs for your comment on the multi-year La Niña simulation capability of the GC-\nLa Niña events. These diagrams show that both the GC-OGCM and ROMS can\nbasically simulate the two peaks of SST cold anomalies; however, the centers of the\nsimulated amplitudes are relatively westward—this is consistent with the common bias\nin climate models where the simulated ENSO amplitude center is shifted westward.\nNevertheless, overall, both ROMS and the GC-OGCM are capable of simulating multi-\nyear La Niña events.\nFigure R3| Time–longitude composite section of anomalous SST averaged within 2°S-\n2°N for the multi-year La Niña events, derived from (a) observation from OBS\n(ERSSTv5), and (b) OGCM simulation, (c) ROMS during 1980-2022. Stippling\nrepresents values above the 90% confidence level by a two-tailed t-test.\nSome texts are added as follows (Lines 235-239; 259-265; 327-333; 410-413):\n“Most studies have confirmed that imposing negative rainfall anomalies over the\nwestern tropical Pacific, with an amplitude comparable to observed interannual\nvariability, can lead to a cooling exceeding –0.3°C in the central-western equatorial\nPacific. This magnitude of cooling is consistent with the SST changes simulated in our\nexperiment (Table S2).”\n“Furthermore, additional experiments with rainfall anomalies applied at different\nlongitudes confirm that the amplitudes of the SST and MLS responses during multiyear\nLa Niña events are generally consistent across different spatial distributions of rainfall\n(Fig. S6). These experiments collectively demonstrate that the effects of negative\nrainfall anomalies from the first and subsequent years superimpose, contributing to the\nmore pronounced cooling in the second-year La Niña.”\n“It is important to note that our findings, while robust within our experimental\nframework (Table S2 and S3), rely on a single OGCM which exhibits certain biases in\nsimulating ENSO amplitude (Fig. S8). In addition, while numerous studies using\ndiverse models have investigated the effects of rainfall on the tropical Pacific, how\nrainfall-induced salinity changes modulate the predictability of multiyear La Niña\nrequires further verification. Future work should employ a broader suite of coupled\nocean-atmosphere models to improve predictability of these prolonged events 52”\n“Both models successfully capture the spatiotemporal evolution characteristics of\nmultiyear La Niña events. (Fig. S8). Comparing these two models enhances the\nrobustness of our analysis results. Detailed information about the specific\nconfigurations of the models and experiments are provided below (Table S3).”\nOther points:\ndata. It makes no sense to include it in such diagnostics.\n=> Thanks for your comments. In this revision, figure related GODAS is removed.\nL123-124: While previous budget analyses did suggest the critical roles of the two\nterms, your Fig. 2a seems not. Can you give some explanations?\n=> We greatly appreciate your correction. We have revised this section. As you can see,\nboth the vertical mixing process and the horizontal advection term play crucial roles in\nthe initial development stage of multi-year La Niña. However, since the vertical mixing\nprocess is a residual term from the budget analysis, we have not overinterpreted this\nterm. Meanwhile, it can be observed that the salinity anomalies induced by the zonal\nmean current in the advection term are highly significant across all three key periods;\ntherefore, we have focused on discussing the contribution of this term. In addition, to\navoid misunderstanding, we have revised the original term \"zonal advection\" to \"ocean\ndynamical processes\" in the Abstract, main text, thereby emphasizing the critical role\nof entire oceanic processes in the initial stage of multiyear La Niña events. Relevant\nrevisions to the budget analysis in Section 2 have also been made to provide clearer\nelaboration (Lines 119-151).\nL167: Other than saying “systematic biases”, it is more informative to specify what\nbiases you are referring to.\n=> Rephrased as follows (Lines in 163-166):\n“However, due to the weaker vertical mixing and entrainment associated with wind\nstirring in the GC-OGCM47, the simulated amplitude of interannual MLS anomalies is\n28\nsubsequently underestimated..”\nL200: Figs 4a and 4b=>4b and 4c;\n=> Corrected as suggested\nL204: winter=>summer; Fig. 4c=>4a\n=> Corrected as suggested\nL209: after “instability”, add “in the western Pacific”\n=> Added as suggested\nL247: Can you explain how you derive the numbers? In Extended Data Table 1, it is\nuseful to include the reference information in parallel to the observed values.\n Thanks for your suggestion, some texts in Table S1 are added as follows:\n“In this study, the Rainfall X1 and Rainfall X2 experiments apply the observed\ncomposite and doubled rainfall negative anomalies, respectively, to investigate their\nimpacts on ocean temperature and salinity. The constant Rainfall experiment imposes a\nsustained anomaly of –9 mm day⁻¹ in the western-central equatorial Pacific (140°E–\n180°, 5°S–5°N) from June(0) to December(1). Observed SST anomalies are derived\nfrom ERSSTv5 (1980–2022). All reported effects are calculated as the differences\nbetween the aforementioned sensitivity experiments and the climatological rainfall\nexperiment under identical composite atmospheric forcing during multiyear La Niña\nevents. Because the Rainfall X2 experiment more realistically captures the observed\namplitude of rainfall-induced salinity variations in the GC-OGCM, we used its results\nto quantify the effect of rainfall on multiyear La Niña amplitude. This amplitude is\n29\ncharacterized by the Niño3.4 index during the D(0)JF(1) and D(1)JF(2) periods. This\napproach partially corrects for the model's inherent underestimation of salinity anomaly\nstrength and its contribution to ENSO evolution.”\nL414: Are your budget analyses based on monthly or daily outputs?\n=> “In this study, we employed an MLS budget derived from monthly mean of ORAS5\nreanalysis to assess the relative contributions of advection, vertical mixing, and\nfreshwater flux (E-P) to the variations in salinity within the mixed layer.”\nPlease explain how the freshwater flux (evaporation minus precipitation) is represented\nin your OGCMs? By natural boundary condition (Huang 1993)?\n Thanks for your suggestion. Some texts are rephrased as follows (Lines 429-432):\n“In the GC-OGCM, freshwater flux entering the ocean influences oceanic processes\nthrough two primary pathways: first, as a virtual salt flux; and second, by modifying\nthe vertical velocity at the top layer via the natural boundary condition, which\nsubsequently modulates the exchange of heat and momentum 26,75..”\nL484/Sect. d: Ensemble information is needed for the experiments.\n=> Thanks for your suggestion. Table S3 is added to describe the model experiments in\nthe main text. As for L484 in the original version, some texts are added in Lines 454-\n459 as follows:\n“Given that the contribution of heat flux to MLS variability is relatively small over\nshort timescales and acts primarily through its modulation of buoyancy flux and density\n(which subsequently influence mixing and circulation), the direct driver of MLS\n30\nvariability in the rainclim experiment is ocean dynamical processes induced by wind\nforcing. Therefore, we did not explicitly distinguish the contributions of heat and\nmomentum fluxes in this experimental framework\nTable S3 Summary of Main Experiments design\nExperiments I: Impacts of Rainfall on mixed layer salinity during multiyear La Nina\nExpeirment Ocean Atmosphere Experiment Experiment\nName Model Forcing configuration Goals\nThe daily\natmospheric Rainfall is prescribed as an\nGC-\nrain forcing data from interannually varying field during\ninter\nOGCM/ROMS Isolating the\nNCEP/NCAR 1979-2022\ncontributions of\nutilized in this\nrainfall to mixed\nstudy are then\nlayer salinity during\napplied to calculate\nGC- Rainfall is prescribed as its multiyear La Nina\nrain the flux field\nclim\nOGCM/ROMS climatological field\nthrough the bulk\nflux formula.\nExperiments II: Impacts of interannual rainfall anomalies on multiyear La Niña\nExpeirment Ocean Atmosphere Experiment Experiment\nName Model Forcing configuration Goals\nBaseline control\nThe observed monthly climatological\nclimatological experiment, which\nrainfall from 1980 to 2022 is\nrainfall GC-OGCM serves to compare\nprescribed in ocean physical model\nexperiments with subsequent\nsimulations\nexperiments.\nThe rainfall\nComposite surface\nThe interannual anomalies of the anomaly was\nwinds from\ncomposited rainfall during multiyear doubled to match\nNCEP/NCAR\nRainfall X2 GC-OGCM La Niña events are multiplied by a the observed MLS\nduring multiyear\nfactor of two and superimposed on changes as well as\nLa Niña events are\nthe climatological rainfall field their corresponding\nused to force the\nimpacts on SST.\nGC-OGCM model\nObserved rainfall negative anomaly\nfor a duration of\nis utilized force the GC-OGCM\ntwo years.\nmodel for a duration of two years For comparison\nRainfall X1 GC-OGCM\nunder composite surface winds from with Rainfall X2\nNCEP/NCAR during multiyear La\nNiña events\nRainfall (-9mm Constant rainfall anomaly (- The selection of this\nGC-OGCM\nday-1) 9mm/day) is imposed into the region and the\n31\nwestern-central equatorial Pacific magnitude of the\n(140°E-180°) during June (0)- rainfall anomaly are\nOctober (1) to solely affect salinity designed to closely\nmatch the observed\nMLS.\nFig1: Which datasets are used for the analyses?\n Added as suggested (Lines 725-725).\n“Datasets: SST from ERSSTv5, MLS from IAP, and rainfall from GPCP.”\nReferences:\n1. Wu, L., Sun, Y., Zhang, J., Zhang, L. & Minobe, S. Coupled ocean-atmosphere\nresponse to idealized freshwater forcing over the western tropical pacific. J. Clim. 23,\n1945–1954 (2010).\n2. Vialard, J. et al. A Model Study of Oceanic Mechanisms Affecting Equatorial\nPacific Sea Surface Temperature during the 1997–98 El Niño. J. Phys. Oceanogr. 31,\n1649–1675 (2002).\n3. Maes, C., Picaut, J. & Belamari, S. Salinity barrier layer and onset of El Niño in a\nPacific coupled model. Geophys. Res. Lett. 29, 59-1-59–4 (2002).\n4. Zhu, J. et al. Salinity anomaly as a trigger for ENSO events. Sci. Rep. 4, 6821\n(2014).\n5. Kang, X., Zhang, R. H. & Wang, G. Effects of different freshwater flux\n32\nrepresentations in an ocean general circulation model of the tropical Pacific. Sci. Bull.\n62, 345–351 (2017).\n6. Yang, S., Lau, K. M. & Schopf, P. S. Sensitivity of the tropical Pacific Ocean to\nprecipitation-induced freshwater flux. Clim. Dyn. 15, 737–750 (1999).\n7. C. Guan, F. Tian, M. J. McPhaden, S. Hu, F. Wang, Zonal Structure of Tropical\nPacific Surface Salinity Anomalies Affects the Eastern and Central Pacific El Niños\nDifferently. Geophys. Res. Lett. 50, 1–10 (2023).\n8. C. Guan, et al., Zonal Structure of Tropical Pacific Surface Salinity Anomalies\nAffects ENSO Intensity and Asymmetry. Geophys. Res. Lett. 49 (2022).\n9. R. Seager, M. B. Blumenthal, Y. Kushnir, An Advective Atmospheric Mixed Layer\nModel for Ocean Modeling Purposes: Global Simulation of Surface Heat Fluxes. J.\nClim. 8, 1951–1964 (1995).\n10. R. Murtugudde, R. Seager, A. Busalacchi, Simulation of the Tropical Oceans\nwith an Ocean GCM Coupled to an Atmospheric Mixed-Layer Model. J. Clim. 9, 1795–\n1815 (1996).\n33"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the authors' great efforts, and they clearly addressed my concerns. Thus, I have no more questions and recommend publication in Nature Communications.", "reviewer_2": "The authors have well addressed my concerns. The paper can be accepted now."}, "author_response": ""}]

Decision letter