Introduction

ENSO prediction skill has not shown obvious improvement over the past few decades, and in fact, there has been a decline at the turn of the 21st century1. The decreased skill is largely attributed to the increased occurrence of CP El Niños2,3,4. Previous theories of ENSO primarily focused on the positive Bjerknes feedback5 and several negative feedback6,7,8,9,10,11. These theories highlight the wind-driven redistribution mechanism of the warm waters, either between the eastern and western equatorial Pacific (known as the zonal advection mechanism), or between the equatorial and off-equatorial regions (known as the thermocline mechanism). These mechanisms are all essentially adiabatic from the oceanic perspective.

However, insights from both numerical models and diagnostics of the parameterized vertical diffusion/entrainment of temperature suggest that diabatic processes around the mixed layer (ML) base (MLB) may play an important role in the development of both CP and EP ENSO12,13,14,15,16,17. Direct turbulence observations demonstrate that entrainment and turbulent mixing just beneath the MLB frequently occur from the western18 to central and eastern equatorial Pacific19. Recently, using approximately 10 years of nearly continuous turbulence observations by moored sensors at 140°W, 0° (Fig. 1a), Warner and Moum20 found positive feedback of subsurface turbulent mixing at this site to the ENSO-related interannual sea surface temperature (SST) changes (represented by the Niño 3.4 region-based Oceanic Niño Index). While representing a great advance in data-based diagnostic of turbulent mixing’s effect on ENSO, the employed observation site is unfortunately not located in the core area of either CP ENSO or EP ENSO (represented by Niño 4 and Niño 3 regions, respectively). Extensive sampling that encompasses at least the central to eastern Pacific equatorial band are required, therefore, to validate this positive feedback mechanism, and to distinguish its effect on EP from that on CP SST changes. Yet, achieving such comprehensive sampling of subsurface turbulent heat flux remains challenging due to practical difficulties in measuring turbulence. In this work, we present an approach to address this problem; particularly, we reveal that the subsurface turbulent mixing can enhances the CP ENSO as positive feedback.

Fig. 1: Data availability of surface and subsurface heat fluxes respectively into and out of the mixed layer in the equatorial Pacific.
figure 1

a Schematic of the inward surface net heat flux to the mixed layer and outward subsurface turbulent heat flux of the mixed layer base (MLB); the shading displays the mean monthly temperature anomalies over the growing months (Sep., Oct. and Nov., SON) of all El Niño events during 2000–2022; the black thick line shows the mean depth of MLB of the equatorial band, calculated from Argo profiles. The regions of Niño 3, Niño 4 and Niño 3.4 are marked. b Argo profile numbers, counted over 5°(zonal)×3°(meridional) boxes over 2005–2022 (left), and the zonal mean Argo profile numbers averaged over 170°W–110°W at 1° meridional grid (right). c The estimated mean inward surface net heat fluxes (penetrative shortwave heat fluxes considered) and outward subsurface turbulent heat fluxes, \({J}_{q}^{s}\) and \({J}_{q}^{{EL}}\), respectively, averaged over 5° (zonal)×6° (3°S–3°N) boxes along the equatorial band (see the text); shading with dashed bounds denotes the corresponding standard deviation calculated over each box during 2005–2022, shading with solid bonds denotes the 95% bootstrap confidence intervals, while error bars denote the standard deviation of monthly anomalies (i.e., the magnitude of interannual variations).

Full size image

Results

The entrainment layer mixing on the equator

We first focus on identification of mixing events within a layer surrounding the MLB. The reason why we focus on such a layer is because it is just beneath the MLB that turbulent mixing directly transfers the heat in the ML downward and hence influences the ML temperature. This layer is referred to as the entrainment layer (EL) henceforth (see Methods for detailed definitions). The mixing events are identified as density overturns21,22,23, because the latter is one of the key features during shear instability-associated billows and turbulences in the shear flow21. The mixing events are determined across the low latitude Pacific Ocean from the long-term (2000–2022), wide-coverage, high-resolution and quality-controlled Argo profile dataset; there, the Argo profiles are concentrated to the equatorial band (Fig. 1b), with more than 40,000 Argo profiles available in the central-to-eastern equatorial Pacific Ocean (170°E–100°W, 3°S–3°N), providing abundant data sources. At each Argo profile, a mixing event is determined only when the density overturn patch covers three or more consecutive vertical grid points, and an EL mixing event is determined only when part or full of a density overturn is lying within the outmost bounds of the EL (see Methods). Combined with a velocity shear-based, wind-driven turbulent mixing scaling scheme as described below, the Argo profiles in this region lead to abundant subsurface turbulent heat flux estimates; further, associated with the widely available surface net heat fluxes (Fig. 1a, c), the diabatic effect on the ML temperature change can be obtained eventually.

Before giving the details of the turbulent mixing scaling scheme, we first show the occurrence of the EL mixing events and analyze their formation mechanisms. The detected density overturns reveal widespread and extensive emergence of EL mixing events throughout the equatorial Pacific Ocean (Fig. 2a), with a concentration of occurrence rate (\({r}_{{ELM}}\)) in the central-to-eastern equatorial Pacific (170°E–110°W, 3°S–3°N), a maximum \({r}_{{ELM}}\) of >50% centered at 160°W (Fig. 2b). The mean \({r}_{{ELM}}\) over this concentration region (170°E–110°W, 3°S–3°N) is 41.2%; here, \({r}_{{ELM}}\) is meridionally consistent from 3°S to 3°S (Fig. 2c). In contrast, \({r}_{{ELM}}\) diminishes to <20% west of 170°E and east of 110°W. Between 140°–170°E (the Pacific warm pool), \({r}_{{ELM}}\) is also meridionally consistent; whereas, it is smaller at the equator than off-equator between 110°–90°W. This distribution pattern suggests that the central-to-eastern equatorial Pacific (170°E–110°W, 3°S–3°N) experiences more frequent diapycnal communications between the surface and subsurface layers than both the far eastern Pacific (east of 110°W) and the western Pacific (WP, west of 170°E) regions.

Fig. 2: The entrainment layer (EL) mixing in the equatorial Pacific.
figure 2

a Occurrence of EL mixing (\({r}_{{ELM}}\)) at 5° × 3° boxes over 2005–2022. b \({r}_{{ELm}}\) along the equatorial band, calculated over 5° (zonal) × 6° (3°S–3°N) boxes; the ranges of Niño 3, Niño 4, and Niño 3.4 are marked. c Zonal mean \({r}_{{ELM}}\) averaged over the western (140°–170°E), central-to-eastern (170°E–110°W) and far eastern (110°–90°W) at 1° meridional grid. d Monthly EL mixing occurence \({r}_{{ELM}}\) and zonal wind stress \({\tau }_{x}\) averaged over the 170°E–110°W and 3°S–3°N region; inlets show the correlation coefficient between the 13-month running smoothed \({r}_{{ELM}}\) and \({\tau }_{x}\). e The parameterized turbulent dissipation rate \({\varepsilon }_{{ELM}}\) versus the observed dissipation rate \({\varepsilon }_{{obs}}\). \({\varepsilon }_{{obs}}\), N (the buoyancy frequency) and S (the vertical shear of zonal velocity) are calculated from daily TIWE data averaged over 20–30 m (see Methods). f The mean \({k}_{{EL}}\) (in the 10-based logarithm scale) averaged over 5° (zonal)×6° (3°S–3°N) boxes along the equatorial band; shading denotes the 95% bootstrap confidence intervals.

Full size image

Firstly, wind stress forcing offers one of the basic interpretations for occurrence of EL mixing. Stronger zonal wind stress (stronger trade wind) leads to higher occurrence of EL mixing via two mechanisms. The first is it provides stronger vertical shear and hence more shear instability condition in subsurface layer, which will be demonstrated in the following subsection, and the second is it leads to stronger and quicker enhancement of entrainment turbulence around the MLB, which are demonstrated by laboratory experiments24, direct numerical simulation25, and observed by several years of turbulence measurements19. Our results show large \({r}_{{ELM}}\) and obvious response of \({r}_{{ELM}}\) to wind stress at multiple time scales in the central-to-eastern equatorial Pacific, because, in this region, both the mean magnitude and variations of the trade wind are strong. In contrast, the WP ( < 160°E) is characterized by near-zero winds and hence few mixing events (Supplementary Fig. 1h). In the southeast of the equator, the mean zonal wind stress is larger than northeast of the equator, which is associated with higher \({r}_{{ELM}}\) there. The wind’s effect on \({r}_{{ELM}}\) is particularly obvious in time series. In the high \({r}_{{ELM}}\) region (170°E–110°W, 3°S–3°N), \({r}_{{ELM}}\) exhibits distinct seasonal variation and ENSO-related interannual variation features; the monthly \({r}_{{ELM}}\) in this region coincides with the magnitude of zonal wind stress (\({\tau }_{x}\)) (Fig. 2d) at these temporal scales. To emphasize the ENSO-focused interannual variations, 13-month running means of both the monthly \({r}_{{ELM}}\) and \({\tau }_{x}\) are calculated, which show a statistically significant correlation (with a peak correlation coefficient R = –0.75, 95% CI = [–0.82 –0.64], and P-value < 0.001 over 2010–2022).

Secondly, local background shear stability condition26 offers another basic interpretation for occurrence of EL mixing. The occurrence of EL mixing reflects fulfillment of local shear instability therein. On average, the subsurface of the central equatorial Pacific (170°E–140°W) features as moderate background vertical shear of zonal velocity from the South Equatorial Current (SEC) and the Equatorial Undercurrent (EUC), and weak stratification due to the Pacific warm pool, while the near-eastern Pacific (140°–110°W) is featured as strong shear and moderated stratification due to the shoaling of the thermocline (Supplementary Fig. 1a–f); as a result, this central-to-eastern equatorial Pacific region (170°E–110°W) has more likelihood to form a state of marginal instability (MI, i.e., the Richardson number Ri ≈ 0.25; Supplementary Fig. 1g), and results in frequent intermittent mixing events27,28. The center of low Ri locates at 160°W, the same as the center of highest \({r}_{{ELM}}\). In contrast, in the far eastern Pacific (east of 110°W), though the shear of zonal velocity is strong, the thermocline is very sharp and hence the stratification is strong, while in the western Pacific (west of 170°E), though the stratification is mostly moderate, the shear is too weak; as a result, the likelihood of shear instability in both regions is greatly reduced, so as the mixing events. Therefore, the zonal pattern of subsurface Ri (Supplementary Fig. 1g) mimics that of \({r}_{{ELM}}\) (Fig. 2b).

It is noted that the roles of wind stress and shear in promoting the generation of EL mixing are not independent, particularly under MI condition which is frequently fulfilled in the central-to-eastern equatorial Pacific region (170°E–110°W). On the one hand, wind stress controls the strength of subsurface background shear and the initiation of entrainment at the MLB, on the other hand, strong background shear or MI condition favors shear instability and the wind-driven entrainment to develop into deep cycle turbulence. Overall, the MI in the subsurface layer provides the precondition for the occurrence of turbulence, while the wind stress initiates the turbulence there.

Accordingly, we develop a method to quantify the mixing coefficients, including turbulent kinetic energy dissipation rate, turbulent diffusivity and thus turbulent heat fluxes, based on both wind stress and local shear. This enables us to reveal the spatial-temporally varying EL mixing across both the central and eastern Pacific regions. We start from a recent work by Whitt et al.29, which proposed a scaling model (hereafter Whitt et al.’s scaling model) for the maximum vertical turbulent buoyancy fluxes (\({F}_{b}\)) of strong mixing in the subsurface layer, expressed as \({F}_{b}\approx 0.2{u}_{*}^{2}S\), where \(S\) is the vertical shear of the horizontal velocity, and \({u}_{*}^{2}\) represents the wind stress-induced friction velocity squared (\({u}_{*}^{2}=|{\tau }_{x}|{\rho }^{-1}\)), with \(\rho\) denoting water density. Appropriately, the detected EL mixing events (overturns) in the present study can be regarded as the ‘strong mixing’ near the MLB, which induces maximum turbulence buoyancy/heat fluxes therein. Furthermore, employing the turbulent kinetic energy equation, and assuming that the turbulence is undergoing a nearly steady state (Ri≈0.25, i.e., \(S\approx 2{{\rm{N}}}\)) and the turbulent kinetic energy is on nearly equilibrium30, we extend Whitt et al.’s scaling model and propose estimation models for both the turbulent kinetic energy dissipation rate (\({\varepsilon }_{{ELM}}\)) and diapycnal diffusivity (\({k}_{{ELM}}\)) of the EL mixing in the equator (within ±3°, 160°E–110°W), i.e.,

$${\varepsilon }_{{ELM}}\approx 1.6{u}_{*}^{2}N,$$
(1)

and

$${k}_{{ELM}}\approx 0.32{u}_{*}^{2}/N,$$
(2)

where \(N\) represents the mean buoyance frequency over the EL. Surprisingly, the schemes show good consistency with observational data at 140°W, 0°N (Fig. 2e; refer to Methods for detailed derivations of Eqs. 1 and 2, as well as uncertainty and validation analysis for these estimation schemes). The advantage of these schemes lies in that, their determinative factors \({u}_{*}^{2}\) and \(N\), as well as the strong mixing condition, can be directly obtained from available wind stress and Argo profile datasets.

Conversely, for weak or non-EL mixing conditions where Whitt et al.’s scaling model is inapplicable, a widely-used background diffusivity \({k}_{{ELnomix}}={10}^{-5}\,{{{\rm{m}}}}^{2}\,{{{\rm{s}}}}^{-1}\) is prescribed (Methods).

Finally, we present the EL mixing diffusivity \({k}_{{EL}}\) for all the Argo profiles in the equatorial region by combing the two conditions:

$${k}_{EL}=\left\{\begin{array}{cc}0.32{u}_{*}^{2}/N,\hfill & strong\,EL\,mixing\,exists \hfill \\ {10}^{-5} \, ({m}^{2}{s}^{-1}),\hfill & no\,EL\,mixing\,detected\end{array}\right.$$
(3)

The obtained mean \({k}_{{EL}}\) (including both \({k}_{{{\rm{ELM}}}}\) and \({k}_{{{\rm{ELnomix}}}}\)) along the equatorial band is shown in Fig. 2f. It exhibits an average magnitude of 10−4 m2 s−1, yet is several times higher in the very central equatorial Pacific region (180°–140°W) compared to either west or east of it.

The determination of \({k}_{{EL}}\) provides the basis for direct estimation of turbulent heat flux in the EL, which is calculated as:

$${J}_{q}^{{EL}}=\rho {C}_{p}{k}_{{EL}}{T}_{z},$$
(4)

where \(\rho\) is water density, \({C}_{p}\) the heat capacity, and \({T}_{z}\) the bulk mean vertical gradient of sorted temperature averaged either over the EL mixing patches for the strong EL mixing category, or over the EL for the non-EL mixing category (Methods; Supplementary Fig. 2a). Here, \({J}_{q}^{{EL}}\) represents the diapycnal turbulent heat flux downward from the mean depth of the EL, \({h}_{{ELM}}\), which closely approximates the depth of the MLB (\({h}_{{MLB}}\)) based on statistics (Supplementary Fig. 2b).

Meanwhile, at each Argo profile ___location, the net inward diapycnal heat flux from the sea surface into the ML, \({J}_{q}^{s}\), is calculated by subtracting a portion of shortwave radiation that penetrates further downward at \({h}_{{ELM}}\) (\({I}^{{h}_{{ELM}}}\)) from the net surface heat flux (\({J}_{q}^{0}\)) (Methods).

As a result, long-term, wide-coverage data of \({J}_{q}^{{EL}}\) and \({J}_{q}^{s}\) are quantified (Supplementary Fig. 3a, b), which is crucial for assessing the diabatic contribution to temperature variations in the ML across various time scales and regions. The temporal means of \({J}_{q}^{{EL}}\) and \({J}_{q}^{s}\) along the equatorial band (within ±3°) show comparable magnitudes, ranging from ~10 Wm−2 at 170°E to ~150 Wm−2 at 110°W, and prominent seasonal and interannual variations (also see Fig. 1c). Both \({J}_{q}^{{EL}}\) and \({J}_{q}^{s}\) exhibit magnitudes and seasonal variations quite like those observed at 140°W, 0°N (Supplementary Fig. 4a; ref. 31). The magnitudes of ENSO-focused interannual variations for both \({J}_{q}^{{EL}}\) and \({J}_{q}^{s}\) reach 50 Wm−2 (Fig. 1c; Supplementary Fig. 3c, d), indicating that \({J}_{q}^{{EL}}\) can exert a great influence on ML temperature comparable to that of \({J}_{q}^{s}\) at the ENSO scale.

The entrainment layer mixing’s response to ENSO

Horizontally, the monthly anomalies of \({J}_{q}^{s}\) (i.e., \({J}_{q}^{{s}^{\prime}}\)) show a clear out-of-phase relationship with the ENSO-related interannual varying SST anomalies (SSTAs), and demonstrate nearly zonal coherence across the entire central-to-eastern Pacific basin (170°E–110°W) (Supplementary Fig. 3c). In contrast, the monthly anomalies of \({J}_{q}^{{EL}}\) (i.e., \({J}_{q}^{{EL}^{\prime} }\)) display an approximately opposite coherent pattern between the CP region (which is defined specifically for 170°E– ~ 140°W and 3°S–3°N) and the EP region (which is defined specifically for ~130°–110°W and 3°S–3°N), although with their border migrating east- and westward year-by-year (Supplementary Fig. 3d).

Note that both the CP and EP regions defined in the present study are different from the conventional Niño 4 and Niño 3 regions (see Fig. 2b, f), respectively. While the Niño 4 and Niño 3 regions represent the core regions on SST anomaly of CP ENSO and EP ENSO, respectively, the defined CP and EP regions here signify the zonal discrepancy of the diabatic effect on the interannual SST changes. We will show that this pattern is associated with the zonal discrepancy of the ENSO-focused interannual variations of wind stresses in the following sections. Zonally, the CP region covers most of the Niño 4 region and also the western part of the Niño 3 region, so that it can represent the conventional CP ENSO region, and can represent the western part of the conventional EP ENSO region as well; in the meanwhile, the EP region as defined in the present study can only represent the eastern part of the conventional EP ENSO region.

In the CP region, the easterly trade winds undergo strong intensification during the growing months (September to November, SON) of the La Niña phase (Fig. 3b, c; Supplementary Fig. 5a), which impart greater momentum to the upper ocean. Moreover, the enhanced easterly winds strengthen the westward SEC in the upper layer, and induces enhanced eastward EUC in the subsurface layer by resulting in west-east pressure gradient32. Consequently, the squared shear (\({S}_{u}^{2}={(\frac{\partial u}{\partial z})}^{2}\)) of the zonal velocity (\(u\)) is amplified within the EL (and further below); meanwhile, the stratification (\({N}^{2}\)) remains almost unchanged in the CP region, or even slightly decreased particularly west of 160°W. Consequently, as displayed in Supplementary Fig. 5c–f, the reduced shear squared, \({S}_{{uRed}}^{2}={S}_{u}^{2}-{4N}^{2}\) (or the reciprocal of the Richardson number, \({{Ri}}^{-1}=\frac{{S}_{u}^{2}}{{N}^{2}}\)), becomes anomalously positive, which provides more favorable shear instability conditions and hence mixing events (Fig. 3e, f). The wind anomalies reverse during the growing months of El Niño phase, which ultimately suppress shear instability and turbulent mixing therein (Fig. 3a, d, c, f; Supplementary Fig. 5a, b).

Fig. 3: Entrainment layer (EL) mixing’s response to ENSO differs between central Pacific (CP) and eastern Pacific (EP) regions.
figure 3

a Monthly zonal wind stress \({\tau }_{x}\) anomalies (vectors) and SST anomalies (shading) (relative to the monthly climatology over 2000–2022) averaged over the growing months (SON) of El Niño years (2002, 2004, 2006, 2009, 2014, 2015, 2018 and 2019). b The same as a, but for La Niña years (2005, 2007, 2010, 2011, 2016, 2017, 2020, 2021, 2022). c The mean \({\tau }_{x}\) of SON along the equatorial band (3°S–3°N) for the aforementioned El Niño, La Niña and all years, respectively; the zonal ranges of CP (170°E–140°W) and EP (130°–110°W) as referred in the present study are marked. d Occurrence of EL mixing (\({r}_{{ELM}}\)) at 5° × 3° boxes for SON of El Niño years (2006, 2009, 2014, 2015, 2018 and 2019). e The same as d, but for SON of La Niña years (2007, 2010, 2011, 2016, 2017, 2020, 2021, 2022). f The same as c, but for \({r}_{{ELM}}\).

Full size image

The analyses above are represented by the differences in \(u\), \({{{\rm{S}}}}_{u}^{2}\),\(\,{N}^{2}\) and \({S}_{{uRed}}^{2}\) between westward and eastward anomalous wind stresses at 170°W, 0°N (Supplementary Fig. 5c–f), where data is available thanks to the TAO project (see Methods). This result also can be inferred from 140°W, 0°N20.

Consequently, more EL mixing is detected during the growing months of La Niña than that of El Niño in the CP region (Fig. 3f; Supplementary Fig. 5b). Therefore, the \({r}_{{ELM}}\) (Fig. 4b), the mean \({k}_{{EL}}\) (Fig. 4c) and the mean \({J}_{q}^{{EL}}\) (Fig. 4d) therein all show strong correlations with the ENSO-related wind stress \({\tau }_{x}\). [In addition, time series also shows the strong correlation between magnitudes of \({\tau }_{x}\) and \({J}_{q}^{{EL}}\) at seasonal time scales (Fig. 4a).] These results clearly illustrate the projection of ENSO-driven variability across the air-sea interface onto the subsurface turbulent mixing.

Fig. 4: Turbulent heat flux (\({J}_{q}^{EL}\)) at the mean depth of entrainment layer (EL) and EL mixing parameters as dependency of zonal wind stress τx in the central Pacific (CP) and eastern Pacific (EP) regions.
figure 4

a Time series of monthly \({\tau }_{x}\) and monthly \({J}_{q}^{{EL}}\) averaged over CP region (170°E–140°W and 3°S–3°N); shading denotes 95% bootstrap confidence intervals. b–d The responses of \({r}_{{ELM}}\), \({k}_{{EL}}\) and \({J}_{q}^{{EL}}\), respectively, to \({\tau }_{x}\); colors show the mean SST monthly anomalies; all properties are averaged over growing months SON; inlets show the correlation coefficients. e–h The same as a–d respectively, but for the EP region (130°–110°W and 3°S–3°N); compared to a–d note the smaller range of \({\tau }_{x}\), as well as the reversed relation between SST monthly anomalies and \({\tau }_{x}\).

Full size image

Compared to the CP region, the EP region exhibits differences in both wind patterns and shear instability conditions during the ENSO cycle. Firstly, the interannual variation of zonal wind stress is weak therein (Figs. 3a–c, 4e), with slight easterly (westerly) anomalies during the growing months (SON) of El Niño (La Niña) phase (Fig. 3a–c), which are opposite to that of the CP region. The inflection point between westward and eastward zonal wind stress anomalies locates somehow at ~135°W. The wind-driven shear squared \({S}_{u}^{2}\) in the EP region experiences a slight increase (decrease) due to the weak westward (eastward) wind anomalies during the growing months (SON) of El Niño (La Niña), following a similar mechanism as in the CP region. It favors (damps) the shear instability from the perspective of shear component. However, the stratification (\({N}^{2}\)) in the EP region becomes weaker (stronger) due to deepening (shoaling) of the thermocline, resulting in slightly favorable (unfavorable) shear instability conditions from the perspective of stratification. As a result, compared to the CP region, there is an obvious but weak increase (decrease) in EL mixing (including \({r}_{{ELM}}\) and the mean \({k}_{{EL}}\); Fig. 4f, g) in the EP region during the growing months of El Niño (La Niña), which reflects the primary driving mechanism of wind stress (Fig. 4e). The inflection point that separates the positive and negative \({r}_{{ELM}}\) interannual anomalies along the equator locates also around 135°W (Fig. 3f), close to the inflection point for the positive and negative wind stress’s interannual anomalies (Fig. 3c). The inflection point determines the border between CP and EP regions as defined in the present study. However, likely due to the reduced vertical gradient, the interannual variation of the mean \({J}_{q}^{{EL}}\) is discrete, thus it shows a less significant relationship with \({\tau }_{x}\) in this region (Fig. 4h).

The entrainment layer mixing’s feedback on ENSO

Now, we examine the diabatic effect on ENSO. A net downward \({J}_{q}^{s}\) warms the ML water, while a net downward \({J}_{q}^{{EL}}\) cools it. Jointly, the difference between \({J}_{q}^{s}\) and \({J}_{q}^{{EL}}\), i.e.,

$$\Delta {J}_{q}={J}_{q}^{s}-{J}_{q}^{{EL}},$$
(5)

represents the overall retained heat via diapycnal heat fluxes in the ML, with the diabatic warming/cooling rate being estimated as20,31:

$${T}_{t} \left|{\atop{\scriptstyle{J}_{q}}}=\frac{\partial T}{\partial t}\right|{\atop{\scriptstyle{J}_{q}}}=\frac{\Delta {J}_{q}}{{\rho }_{0}{C}_{p}{h}_{{ELM}}}.$$
(6)

Similarly, a positive (negative) monthly anomalous \({{J}_{q}^{s}}^{\prime}\) results in anomalous diabatic warming (cooling) of the ML water, while a positive (negative) monthly anomalous \({{J}_{q}^{{EL}}}^{\prime}\) will lead to anomalous diabatic cooling (warming). Meanwhile, their differences, denoted as \(\Delta {{J}_{q}}^{\prime}\) and \({T}_{t}{|}_{{{J}_{q}}^{\prime}}=\frac{\Delta {{J}_{q}}^{\prime} }{{\rho }_{0}{C}_{p}{h}_{{ELM}}}\), respectively, represent the overall anomalous diabatic warming/cooling effects on the ML water.

In addition to the outstanding annual cycle of \(\Delta {J}_{q}\), which governs the seasonal variation of the SST in both the CP and EP regions (Fig. 5a; Supplementary Fig. 4a–c; ref. 31), the interannual variation in \(\Delta {{J}_{q}}^{\prime}\) also stands out (Fig. 5b; Supplementary Fig. 3c, d), which show potential relation to the ENSO cycles (Fig. 5c)

Fig. 5: Difference between the heat flux retained in the mixed layer from sea surface (\({J}_{q}^{s}\)) and the turbulent heat flux at the mean depth of entrainment layer (\({J}_{q}^{EL}\)), i.e., \(\Delta {J}_{q}(={J}_{q}^{s}-{J}_{q}^{EL})\), and its monthly anomaly, \(\Delta {{J}_{q}}^{\prime}\), as dependency of time and longitude along the equatorial band.
figure 5

a, b monthly mean \(\Delta {J}_{q}\) and its monthly anomalies (\(\Delta {{J}_{q}}^{\prime}\)) averaged over 5° (zonal)×6° (3°S–3°N) boxes along the equatorial band; c, the ocean Niño index (ONI); dashed lines denote –0.5 °C and 0.5 °C.

Full size image

We first explore the qualitative diabatic warming/cooling effect on the growing of El Niño and La Niña in both the CP and EP regions. On average, during the growing months (SON) of El Niño phases (Fig. 6a), both \({{J}_{q}^{{EL}}}^{\prime}\) and \({{J}_{q}^{s}}^{\prime}\) are negative in the CP, however, \({{J}_{q}^{s}}^{\prime}\) is smaller in magnitude compared to \({{J}_{q}^{{EL}}}^{\prime}\), which thus jointly lead to anomalous diabatic warming (\(\Delta {{J}_{q}}^{\prime} \, > \, 0\)) of the ML water. Similarly, during the growing months (SON) of La Niña phases (Fig. 6b), both \({{J}_{q}^{{EL}}}^{\prime}\) and \({{J}_{q}^{s}}^{\prime}\) are positive, but the former outweighs the latter, jointly leading to anomalous diabatic cooling of the ML water (\(\Delta {{J}_{q}}^{\prime} \, < \, 0\)). The results demonstrate the positive role of the diabatic warming/cooling mechanism in promoting the development of both El Niño and La Niña in the CP region. This positive role is primarily induced by \({{J}_{q}^{{EL}}}^{\prime}\).

Fig. 6: Entrainment layer (EL) mixing’s feedback on ENSO differs between central Pacific (CP) and eastern Pacific (EP) regions.
figure 6

a Temporal mean of the monthly anomalies \({{J}_{q}^{{EL^{\prime}}}}\) and \({{J}_{q}^{s\,^{\prime}}} \) averaged over the growing months (SON) and over 5° (zonal) × 6° (3°S–3°N) boxes along the equatorial band of El Niño years (2006, 2009, 2012, 2014, 2015, 2018 and 2019); error bars denote the 95% bootstrap confidence intervals; pink and blue shading denote anomalously diabatic warming and cooling to the ML, respectively. b The same as a but for La Niña years (2007, 2010, 2011, 2016, 2017, 2020, 2021, 2022). c Three-month smoothing change rate of monthly SST anomalies (i.e., \(\frac{\partial {SST}^{\prime}}{\partial t} \), x-axis) averaged over the growing months (SON) of El Niño years (2006, 2009, 2012, 2014, 2015, 2018 and 2019, orange) and La Niña years (2007, 2010, 2011, 2016, 2017, 2020, 2021, and 2022, blue) of the CP region (170°E–140°W, 3°S–3°N), versus the anomalies of diabatic warming/cooling rate [\((\frac{\partial T}{\partial t}{|}_{\Delta {J}_{q}})^{\prime}=\frac{\Delta {{J}_{q}}^{\prime} }{{\rho }_{0}{C}_{p}{h}_{{ELM}}}\))] averaged over the same months and region (filled circles, y-axis), as well as versus the anomalies of diabatic warming/cooling rate caused purely by \({{J}_{q}^{{EL}}}^{\prime}\), i.e., \((\frac{\partial T}{\partial t}{|}_{{J}_{q}^{{EL}}})^{\prime}=\frac{{{J}_{q}^{{EL}}}^{{\prime} }}{{\rho }_{0}{C}_{p}{h}_{{ELM}}}\) (filled squared, y-axis) and versus the anomalies of diabatic warming/cooling rate caused purely by \({{J}_{q}^{s}}^{\prime}\), i.e., \((\frac{\partial T}{\partial t} |_{{J}_{q}^{s}})^{\prime}=\frac{{{J}_{q}^{s}}^{{\prime} }}{{\rho }_{0}{C}_{p}{h}_{{ELM}}}\) (diamond, y-axis); confidence intervals are the 95% bootstrap confidence intervals of the monthly data. Like in ref. 20, only those with monthly SST anomalies over 0.15 °C for El Niño phase and less than –0.15 °C for La Niña phase are considered. d The same as c, but for the EP region (130°–110°W, 3°S–3°N).

Full size image

In the meanwhile, in the EP region, \({{J}_{q}^{s}}^{\prime}\) is significantly negative with a large magnitude during the growing months (SON) of El Niño phase, leading to anomalous cooling of the ML water therein, while \({{J}_{q}^{{EL}}}^{\prime}\) is weak and positive (significant only east of 120°W), also contributing to anomalous warming. Together, they induce significant anomalous diabatic cooling (\(\Delta {{J}_{q}}^{\prime} \, < \, 0\)) (Fig. 6a). In contrast, during the growing months of La Niña phase, the strong positive \({{J}_{q}^{s}}^{\prime}\) and weak \({J}_{q}^{{EL}^{\prime} }\) lead to anomalous warming (\(\Delta {{J}_{q}}^{\prime} \, > \, 0\)) therein (Fig. 6b). It is obvious that this overall diabatic effect plays to repress the growth of both the El Niño and La Niña in this EP region, and is controlled by the \({{J}_{q}^{s}}^{\prime}\) component.

Now let’s quantify their contributions to the SST change. In the CP region (Fig. 6c), during the growing months of El Niño, the mean SSTA increases at a rate of ~+0.2 °C/month (also see Supplementary Fig. 6g), while the diabatic warming mechanism contributes to ~+0.35 °C/month (also see Supplementary Fig. 6e), obviously outweighing other factors and dominating the warming trend. Similarly, during the growing months of La Niña, the mean SSTA changes at a rate of ~–0.14 °C/month), whereas the diabatic cooling mechanism contributes to ~–0.29 °C/month, once again dominating the cooling trend. Notably, these contributions primarily stem from the larger ENSO-related \({J}_{q}^{{EL}^{\prime} }\), rather than \({J}_{q}^{s\,^{\prime} }\) (Fig. 6c; Supplementary Fig. 6a, c,e), which is associated with rates of +0.48 and –0.15 (–0.50 and +0.22) °C/month, respectively, for the development of El Niño (La Niña). Moreover, the positive contribution from the diabatic warming/cooling mechanism exceeds the positive contribution from the horizontal advection in the CP region, both of which should be balanced by the negative contribution from the vertical advection and horizontal diffusion terms (Supplementary Fig. 7a, b).

In contrast, the diabatic warming/cooling mechanism exhibits a negative effect on the development of either the El Niño or La Niña in the EP region (Fig. 4b, c). During the growing months of El Niño, the ML water warms at a rate of ~0.21 °C/month, while \({J}_{q}^{s\,^{\prime} }\) and \({J}_{q}^{{EL}^{\prime} }\) lead to anomalous cooling at a rate of ~–1.0 °C/month and warming at a rate of ~+0.24 °C/month, respectively. Obviously, while \({J}_{q}^{{EL}^{\prime} }\) promotes the growth of El Niño, \({J}_{q}^{s\,^{\prime} }\) plays a dominant damping role; in sum, the diabatic mechanism represses the growth of the El Niño therein. Similarly, during the growing months of La Niña, the ML water cools at a rate of ~–0.21 °C/month, while \({J}_{q}^{{EL}^{\prime} }\) leads to anomalous cooling of ~–0.18 °C/month, \({J}_{q}^{s\,^{\prime} }\) leads to diabatic warming of ~+1.5 °C/month, dominant the damping of the La Niña here (also see Supplementary Fig. 6f, h).

This discrepancy relative to the CP region arises due to that \(\Delta {J}_{q}^{\prime}\) is dominated by \({{J}_{q}^{s}}^{\prime}\) in this region (Fig. 6a, b, d; Supplementary Fig. 6b, d, f), which varies inter-annually out-of-phase with SSTA, while \({J}_{q}^{{EL}^{\prime} }\) varies weakly, discrete along with wind stress (Fig. 4h), and only significant east of 120°W, hindering its diabatic effect on the ENSO development in the EP region. In other words, the diabatic warming/cooling mechanism acts to suppress the development of ENSO in the EP region (Supplementary Fig. 6f); this negative contribution is reversed by the positive contribution from both the horizontal and vertical advection in this region (Supplementary Fig. 7c, d).

In summary, the diabatic warming/cooling mechanism plays different roles in the two regions. In the CP region, it serves as a prominent positive mechanism to the development of both El Niño or La Niña. In turn, this positive mechanism is intrinsically tied to the EL mixing, which is further primarily governed by the ENSO-associated wind stress variations. Whereas, in the EP region, it represses the development of ENSO primarily through the \({J}_{q}^{s\,^{\prime} }\) component.

Discussion

Note that, in the above section, we analyzed the diabatic effect on the EP and CP regions, rather than on the conventional EP and CP ENSOs. This is because we don’t have enough EP ENSO samples during the Argo era. We encountered 7 El Niños since this century (2002/03, 2006/07, 2009/10, 2014/15, 2015/16, 2018/19, 2023/24), but only the later 5 El Niños occurred during the time of enough Argo profiles. Moreover, among the 5 El Niños, none of them belongs to a pure EP type, instead, the 2009/10, 2014/15, and 2018/19 El Niños belong to the CP type, while the 2015/16 and 2023/24 El Niños belong to the mixed-type. Therefore, it is not possible to distinguish CP El Niño from EP El Niño with the present data.

However, the present analysis on CP region is analogues to the CP ENSO. The 2 mixed-type El Niños also showed a CP ingredient, which manifests as cores of positive SST anomaly and eastward zonal wind anomalies at the CP region; in the meanwhile, they also showed a weak westward zonal wind anomaly in the EP region, resembling the 3 CP El Niños. These features indicate that the pattern of the controlling factor for the EL mixing, i.e., the wind stresses, of the El Niños are in general similar to that of a typical CP El Niño. The same happens to the La Niñas: during all EP type, mixed-type or CP type La Niñas happened during the Argo era, the wind stresses always showed a westward zonal anomaly in the CP region and weak (insignificant) eastward zonal anomaly in the EP region, resembling that of a typical CP type La Niña. Consequently, we conclude that, the derived diabatic positive feedback mechanism is applicable to CP ENSO, therefore, we would call it a diabatic CP ENSO positive feedback mechanism. The positive feedback loop is schematized in Fig. 7a, b. Working with the conventional Bjerknes positive feedback, it propels the growth of SST anomaly in the CP region till counteracted by negative feedback mechanisms.

Fig. 7: Schematic of the diabatic central Pacific (CP) ENSO dynamics and the positive feedback to CP ENSO.
figure 7

a, b represent the El Niño and La Niña growing stages, respectively; the contrast diabatic effects on the eastern Pacific (EP) region is also shown. Shading on the x-z section denotes the monthly \({N}^{2}\) anomalies (\({N}^{2{\prime} }\)) averaged over SON of El Niño and La Niña years since 2000, respectively, and contours are corresponding isotherms (in °C), based on IAP data (see Methods).

Full size image

In the meanwhile, since the CP region defined in this study covers the western part of the conversional region of EP ENSO, the diabatic positive feedback is applicable to the western part of the EP ENSO, yet not to its eastern part, so that we can’t call the diabatic effect as a positive feedback mechanism to the EP ENSO – it becomes only a regional effect to the EP ENSO.

This work also reveals the diabatic effect on the decaying of the ENSO phases. In general, during the decaying months (January to March) of either El Niño or La Niña events, the diminishing wind stress anomalies in the CP region weaken the positive feedback effect of \({{J}_{q}^{{EL}}}^{\prime}\), eventually hindering its ability to further amplify the temperature anomalies (Supplementary Fig. 6e). The weakened wind stress anomalies stem from the decaying of the ENSO. The results again suggest that the revealed diabatic positive feedback is adherent to and complementary to the Bjerknes positive feedback, because the key process during both the feedback is the ENSO related wind stresses.

Several studies have revealed the importance of ocean mixing in shaping the temperature variation and pattern of the equatorial Pacific. For example, researchers have shown that subsurface mixing exerts control effect over the seasonal variations of the equatorial Pacific SST31 (Supplementary Fig. 4a), and that thermocline mixing helps mitigate systematic SST biases between model simulations and observational data in the tropical Pacific33 (We contrast the local nature of the impact of EL mixing to the more remote impact of mixing in the thermocline33.) Building upon the pioneer work of Warner and Moum that is resulted from turbulence observations at a single site20, the present study capitalizes on the long-term and widely covered Argo data, again signifying the importance of subsurface mixing to climate; in addition, the present study also reveals the difference in subsurface mixing’s effect on CP and EP ENSOs.

Nowadays, the parameterization of sub-grid diapycnal mixing remains one of the foremost uncertainties in coarse-resolution numerical models. There is optimism that the spatial and temporal variation of the EL mixing could be incorporated into the climate models19,34. Doing so holds the promise of producing ENSOs with appropriate magnitudes and types, thereby enhancing our understanding and predictive capabilities of this pivotal climate phenomenon.

Methods

The Argo profile data

The Argo (Array for Real-time Geostrophic Oceanography) profile data of year 2005–2022 (http://www.argodatamgt.org/Access-to-data/Argo-GDAC-ftp-and-https-servers) is employed. The data provides profiles of temperature, salinity and corresponding pressure above 2000 m in the ocean, with an accuracy of ± 0.002°C for temperature, ± 0.005 psu for salinity, and ± 2 decibar for pressure35. Only those quality-controlled profiles with a vertical resolution of 2 m or less are used. Totally >160,000 Argo qualified profiles are available in the tropical Pacific Ocean (Fig. 1b).

The ERA5 reanalysis data

The zonal wind stress (\({\tau }_{x}\)), net surface heat flux (\({J}_{q}^{0}\)) and short-wave solar radiation (\({I}^{0}\)) components of this dataset are used for the estimation of \({J}_{q}^{{EL}}\) and \({J}_{q}^{s}\) (see below). The ERA5 data is the fifth generation of ECMWF reanalysis36, covering the period of 1940 to present. ERA5 is of high spatial (0.25°) and temporal (hourly) resolutions. An obvious improvement is observed in the ERA5 over previous versions. This dataset is available at https://doi.org/10.24381/cds.adbb2d47. Daily means are used for analysis, and the variables are processed onto the positions and dates of the Argo profiles.

The long-term TAO data

TAO is the abbreviation for the Tropical Atmosphere and Ocean project37. This dataset provides ocean velocity, temperature, salinity profile data and wind stress in the central to eastern tropical Pacific, available at http://www.pmel.noaa.gov/tao/data_deliv. Daily ocean data at 170°W, 0°N and wind stress data in the vicinity are used in determination of conditions of anomalous westward and eastward winds near 170°W, 0°N. The anomalous westward and eastward zonal wind conditions are obtained during 2008–2020 and over 180°–160°W and 5°S–5°N (a 20° × 10° region centered at 170°W, 0°N).

The IAP data

This dataset consists of monthly temperature and salinity data from 1960 to present, with a horizontal resolution of 1° × 1° and 41 vertical levels from the surface to 2000 m38. It is used for calculating the \({N}^{2}\) anomalies in Fig. 7. The data is freely downloaded from: http://www.ocean.iap.ac.cn/.

The GHRSST data

The Level 4 of Group for High Resolution SST (GHRSST; https://www.ghrsst.org/)39 is adopted for calculating the SST change rate during El Niño and La Niña phases in both the CP and EP regions (Fig. 6c, d; Supplementary Fig. 6g, h). The L4 gridded products are generated by combining complementary satellite and in situ observations within Optimal Interpolation systems. This dataset provides daily SST at 0.25° horizontal resolution, and is downloaded from: ftp://ftp-oceans.ncei.noaa.gov/pub/data.nodc/ghrsst/L4/GLOB/JPL/MUR/.

The TIWE data

This data refers to the turbulence measurements during the Tropical Instability Wave Experiment (TIWE) during the fall of 1991 at 140°W, 0°N40. Nearly 6000 casts of microstructure temperature, conductivity, and shear measurements in the upper 200 m were made. Simultaneous horizontal velocity was measured by the ship-mounted Acoustic Doppler Current Profilers (ADCPs). The data was provided by Dr. Lien via personal communication. The data were averaged daily into 36 profiles with the vertical resolution of 1 m from November 4 till December 12.

The Johnson et al. (2002) data

This dataset refers to the potential temperature, salinity, potential density, and velocity data that are obtained in the equatorial band41, which is downloaded from https://floats.pmel.noaa.gov/sites/default/files/atoms/files/meanfit_m.cdf_.zip via Dr. Johnson’s website https://floats.pmel.noaa.gov/gregory-c-johnson-home-page.

The OSCAR data

This dataset refers to the Ocean Surface Current Analyses Real-time (OSCAR), which contain ocean surface mixed layer velocities calculated from satellite-sensed sea surface height gradients, ocean vector winds, and SST fields using geostrophy, Ekman, and thermal wind dynamics42. The v2 is employed, which has a horizontal resolution of 1/3° and temporal resolution of 5 days for the period of 1993–present. It is downloaded from: http://apdrc.soest.hawaii.edu/erddap/griddap/hawaii_soest_bb96_63ae_8950.html.

Definitions of the mixed layer (ML) and entrainment layer (EL) of Argo profile data

The ocean surface ML is a physically well-defined layer and can be determined by many previously proposed criteria. In this study, the ML is determined such that the density at its base is 0.01 kg/m3 greater than the near-surface (at 10 m depth)43. The EL is a layer that covers both the lower flank of the ML and the upper flank of the pycnocline, so that internal waves and shear instability can exist, which thus induce entrainment, turbulence and finally diapycnal heat transport across the ML base (MLB). However, no commonly applicable criteria were available yet for its definition. In this study, for the sake of objectivity, we define a layer of 20 m thick as the EL, with its top being 5 m above and its bottom being 15 m below the MLB (Supplementary Fig. 2a). This definition is achieved also by considering the following reasons. Firstly, the purpose of determining an EL is to detect turbulent mixing that may induce cross-MLB turbulence, so that the vertical boundaries of EL should not be far away from the MLB. Secondly, observations show that, in either the EP29 or western equatorial Pacific18,23, maximum mixing coefficient or turbulent buoyancy/heat fluxes often occur within 30 m below the MLB; while off the equator (at 3°N), large eddy simulation (LES) with as realistic as possible oceanic and atmospheric conditions shows that maximum heat flux often occur within 10 m above the MLB29. Given those observational evidences, and considering otherwise that distant mixing events may not impact the ML, we choose the half distances (5 m and 15 m respectively) to the MLB as the outmost upper and lower bounds of the EL, which is schematized in Supplementary Fig. 2a.

Determination of mixing events in the EL and the mean depth of EL mixing (\({h}_{{ELM}}\))

In the present study, we adopt the density inversion (or density overturn) detection method21 to the long-term (2000–2022), widely covered, high-resolution quality-controlled Argo profiles (see above), to detect mixing events in the EL. An EL mixing event is determined only when the overturn patch covers three or more consecutive vertical grid points, and part or full of it is lying within the outmost bounds of EL (Supplementary Fig. 2a). Since the vertical resolution of the Argo profile data is about 1–2 m, the minimum density inversion scale that can be resolved is 4 m; as such, only the large-scale overturns are detected, which represent ‘strong mixing’ therein. Small scale mixing events possibly exist but can’t be detected. Therefore, it is reasonable to assume that the depths with detected mixing events is associated with largest turbulent fluxes.

The mean depth of EL mixing for each Argo profile, denoted as \({h}_{{ELM}}\), should also be determined. If multiple density overturns were detected within the EL, the mean depth of the largest overturn is defined as the \({h}_{{ELM}}\); if all the density overturns have the same sizes, the mean depth of all of them is defined as the \({h}_{{ELM}}\). Whereas, when no density overturn was detected within the outmost EL bounds, the MLB is defined as the \({h}_{{ELM}}\) (Supplementary Fig. 2a). The mean \({h}_{{ELM}}\) along the equator for the two categories are shown in Supplementary Fig. 2b.

The proposed scaling schemes for estimation of the EL mixing coefficients \(\varepsilon\) and \(k\) based on Whitt et al.’s scaling scheme

Intense EL turbulence is frequently generated by surface wind and buoyancy forcing under sufficient shear conditions therein. Whitt et al.29 proposed a scaling scheme to represent daily maximum turbulent buoyancy and heat fluxes in the subsurface layer in terms of bulk vertical shear and wind stress for strong mixing (whose maximum buoyancy flux roughly is larger than 10−7.5 m2s−3). Theoretically, the turbulent kinetic energy (K) in a shear layer just beneath the MLB evolves due to shear production (\({F}_{m}S\)) and dissipation (\((\varepsilon )\) plus buoyancy flux (\({F}_{b}\)) in the shear layer, while the shear (\(S\)) is driven by the vertical divergence of momentum fluxes, which is written as26,30:

$$\frac{\partial K}{\partial t}={F}_{m}S-\varepsilon -{F}_{b},$$
(7)
$$\frac{\partial S}{\partial t}=\frac{1}{{Hh}}[{u}_{*}^{2}-{F}_{m}\left(h\right)],$$
(8)

where \(h\) is the depth of MLB, \(H\) is the thickness for bulk vertical shear calculation, \(S\) is the bulk vertical shear, \({F}_{m}\) is the vertical momentum flux, and \({u}_{*}^{2}=\frac{|\tau |}{\rho }\) is the friction velocity squared, with \(\tau\) the wind stress and \(\rho\) the water density. Based on theories, observational evidences and large eddy simulation (LES) results, Whitt et al. found that the following relationship,

$${F}_{b}\approx 0.2{u}_{*}^{2}S,$$
(9)

is a plausible scaling model which can explain most of the daily variance of \({F}_{b}\) at the depth of maximum \({F}_{b}\). This model directly connects the subsurface turbulent buoyancy flux and both surface wind forcing and local velocity shear condition.

However, the scaling model (Eq. 9) is difficult to apply to Argo profile data because the latter does not contain information of velocity shear. Therefore, we conduct further derivations to make this model applicable to the abundant Argo profiles. Theory (Eqs. 7 and 8) suggests that when the shear (\(S\)) and turbulent kinetic energy (K) are in a nearly steady state30, which is a good assumption at the daily time scale, then \({u}_{*}^{2}\approx {F}_{m}\left(h\right)\) and \({\varepsilon \approx F}_{m}S-{F}_{b}\); additionally, adopting Eq. 9, we obtain:

$$\varepsilon \approx 0.8{u}_{*}^{2}S.$$
(10)

As such, this scaling model directly relates the turbulent kinetic energy dissipation rate (\(\varepsilon\)) in the subsurface layer to the sea surface wind stress (\(\tau\)) and the subsurface vertical velocity shear (\(S\)). This scaling model is consistent with a theoretical derivation30, and is similar to a previous scaling model that is derived for a thicker layer and hence has a smaller coefficient44.

In addition, note that when the turbulent kinetic energy (K) is in a steady state, the gradient Richardson number (\({Ri}\)) over the shear layer should be around the critical value, i.e.,

$${Ri}=\frac{{N}^{2}}{{S}^{2}}\approx 0.25,$$
(11)

otherwise, K will either increase (when \({{\rm{Ri}}} \, < \, 0.25\)) or decrease (\({{\rm{Ri}}} \, > \, 0.25\))45. A state of \({Ri}\approx 0.25\) is called marginal instability (MI46), which is particularly true in the upper layer of the central equatorial Pacific27,28. Here, \(N=-\frac{g}{{\rho }_{0}}\frac{\partial \rho }{\partial z}\) is the buoyancy frequency and denotes stratification, with \(g\) the gravitational coefficient, \({\rho }_{0}\) = 1025 kgm−3 the reference density, and \(\rho\) the potential density. In a stratified state, a background vertical velocity shear plays as a necessary condition for MI27. Equation 11 thus leads to \(S\approx 2N\), and Eq. 10 is finally transformed to

$$\varepsilon \approx 1.6{u}_{*}^{2}N.$$
(12)

The estimation scheme (Eq. 12) directly relates \(\varepsilon\) in the subsurface layer to the wind stress \(\tau\) and subsurface buoyancy frequency \(N\), which can be more easily obtained from current available datasets.

Furthermore, the diffusivity (\(k\)) can easily be calculated as:

$$k=\varGamma \varepsilon /{N}^{2},$$
(13)

with a widely used constant mixing efficiency, Γ = 0.247. Using scheme Eq. 12, we further rewrite Eq. 13 as:

$$k=0.32{u}_{*}^{2}/N.$$
(14)

This scheme (Eq. 14) builds the relationship between the subsurface diffusivity \(k\) and both \(\tau\) and \(N\).

Note that, according to Whitt et al.’s scaling model, the obtained \(\varepsilon\) and \(k\) should be representative for the position of maximum \({F}_{b}\), while \(N\) and \(S\) are bulk averaged over a certain thickness in the vicinity of the maximum \({F}_{b}\). Therefore, both the depth of maximum \({F}_{b}\) and the thickness for the calculation of the bulk mean need to be further determined.

Firstly, for the case that Argo profiles are with detected EL mixing events, the mixing events are assumed as strong mixing and leads to maximum \({F}_{b}\) therein, because otherwise mixing events won’t be detected, let alone turbulent fluxes. In this case (i.e., Category 1, when strong EL mixing exists, as shown in Supplementary Fig. 2a), the mean depth of EL mixing, i.e. \({h}_{{ELM}}\), is naturally considered as the depth of the maximum \({F}_{b}\), and also of the calculated \(\varepsilon\) and \(k\), again because this position is where the largest overturns are detected. The thickness of EL mixing (\({H}_{{ELM}}\)) is defined as a 10-m thick layer that centered at \({h}_{{ELM}}\); but, if necessary, is further constrained by the outmost bounds of the EL as defined above (Supplementary Fig. 2a). In this category, the diffusivity is calculated as the same as Eq. 14, i.e.,

$${k}_{{ELM}}=0.32{u}_{*}^{2}/N,$$
(15)

where \(N\) is calculated from the sorted density over \({H}_{{ELM}}\).

Accordingly, the dissipation rate can be calculated following the same as Eq. 12, i.e.,

$${\varepsilon }_{{ELM}}=1.6{u}_{*}^{2}N,$$
(16)

and also can be calculated by Eq. 10 when the shear data is available, i.e.,

$${\varepsilon }_{{ELM}}=0.8{u}_{*}^{2}S.$$
(17)

The calculated \({\varepsilon }_{{ELM}}\) with Eqs. 16 and 17 shows consistency and matches the observations at 140°W, 0°N very well (Fig. 1d).

Whereas, it should be noted that, in many times, the wind might be weak, and/or the shear stability condition is not fulfilled for entrainment and turbulent mixing to happen near the MLB, or only small turbulence occurs. In such case, mixing events can’t be detected in the EL on the Argo profiles (Case 2, Supplementary Fig. 2a), the diffusivity is then is assumed to be small and prescribed a widely used background value, i.e.,

$${k}_{{ELnomix}}={10}^{-5}\,{{{\rm{m}}}}^{2}{{{\rm{s}}}}^{-1}.$$
(18)

Here, an EL is still defined (a 10 m thick layer covering the MLB) for the calculation of turbulent heat fluxes.

Finally, the diffusivity of the EL for each Argo profile can be described as:

$${k}_{EL}=\left\{\begin{array}{c}0.32{u}_{\ast }^{2}/N,\qquad \,{{\rm{strong}}}\,{{\rm{EL}}}\,{{\rm{mixing}}}\,{{\rm{exists}}}\,\\ {10}^{-5}\,{{{\rm{m}}}}^{2}{s}^{-1},\,\,\,\,\,\,\,no\,{{\rm{EL}}}\,{{\rm{mixing}}}\,{{\rm{detected}}}\end{array}\right.$$
(19)

In the present study, \({u}_{*}^{2}\) is calculated from the daily ERA5 zonal wind stress. Statistics show that, subject to months and spaces, the occurrence of EL mixing \({\gamma }_{{ELM}}\) ranges between 20% ~ 60%, and the 90 percentile of \({k}_{{ELM}}\) ranges between 10−2.8 ~ 10−2 m2s−1 with the median of 10−2.5 m2s−1 (Supplementary Fig. 8a), which is more than 2 orders larger than \({k}_{{ELnomix}}\). Taking both \({k}_{{ELM}}\) and \({k}_{{ELnomix}}\) into account (i.e., \({k}_{{EL}}\)), the median of them reduces to 10−4.0 m2s−1. \({\varepsilon }_{{ELM}}\) of the EL mixing cases shows a mean magnitude of \({10}^{-8}\) Wkg−1 (ranging from 10−8.5 to 10−7.2 Wkg−1), and is larger in the eastern than central and western equatorial Pacific. Those parameters show clear seasonal and interannual variations (Supplementary Fig. 8b). All those values and features fall into the range of limited observations.

We note that, the schemes Eqs. 17 and 19 are only applied to 3°S–3°N and 170°E–110°W, where the background vertical shear near the MLB commonly exists and hence Whitt et al.’s scheme (Eq. 9)29 and the assumption of MI are valid to a large extent (Supplementary Fig. 1g).

Estimation of the turbulent heat flux (\({J}_{q}^{{EL}}\)) at the mean depth of EL mixing

For each Argo profile, either with or without detected EL mixing events, \({J}_{q}^{{EL}}\) is calculated as \({J}_{q}^{{EL}}=\rho {C}_{p}{k}_{{EL}}{T}_{z}\), where \({C}_{p}\) is the heat capacity, \({T}_{z}\) is the bulk mean temperature gradient averaged over \({H}_{{ELM}}\) for the \({k}_{{ELM}}\) case, and also over a 10-m thick layer around the \({h}_{{ELM}}\) (i.e., the \({h}_{{MLB}}\)) for the \({k}_{{ELnomix}}\) case (Supplementary Fig. 2b). The latter on average is two orders smaller than the former. The calculated \({J}_{q}^{{EL}}\) is prescribed at \({h}_{{ELM}}\). The estimated \({J}_{q}^{{EL}}\) shows a quite reasonable annual cycle (Supplementary Fig. 4a), indicating that the scaling schemes for \({k}_{{EL}}\) and \({J}_{q}^{{EL}}\) are validated.

Uncertainties and verification of the above \({\varepsilon }_{{ELM}}\) and \({J}_{q}^{{EL}}\) estimations

We are aware that several uncertainties exist in the estimation models of both \({\varepsilon }_{{ELM}}\) and \({J}_{q}^{{EL}}\):

  1. (1)

    According to several previous studies as cited above, to avoid extreme fluctuations, the mixing coefficients are usually estimated at daily time scales as done here, however, due to lack of data, the snapshot \(N\) and \({T}_{z}\) obtained from Argo profile are used as the daily proxies. The influence of this assumption is unknown, but it is expected that the adopted temporal and spatial average may smoothed out the bias.

  2. (2)

    The sea surface forcing variables (wind stress and surface heat flux components) are adopted from the ERA5 dataset, which has the longest and finest coverage, shows good performance in the equatorial region48, and represents one of the best available reanalysis data. However, any remaining uncertainty of this dataset in the study region might project to the estimation of the proposed estimation models.

  3. (3)

    The extension of Whitt et al.’s scaling model from 140°W, 0°N and 140°W, 3°N to the broader equatorial band is not direct evidence-based, but by reasonable deduction with indirect information. The Whitt et al.’s scaling model is based on observations at 140°W, 0°N and large eddy simulation (LES) results at both 140°W, 0°N and 140°W, 3°N. The success of the scaling model relies on common existence of vertical shear in the upper layer and marginal instability (MI) condition in the subsurface layer. This condition is likely to be fulfilled alone the narrow equatorial band (3°S–3°N) and in the central-to-eastern region (170°E–110°W) due to common existence of both vertical velocity shear and not-extremely-strong stratification therein, but no direct observational or simulation evidence available to our knowledge. However, a lot of information indicate that MI could be occasionally occur. Based on linear stability analysis, Zhang et al.28 demonstrated that MI occupies a large portion of time in the upper layer in this zonal range at the equator. In the off-equatorial region, the westward flowing SEC, which covers a meridional range of about 5°S to 5°N, provides the background vertical shear in the upper ~50 m41. The SEC is the primary current that provides the shear at the MLB and EL depths, and is enhanced by the EUC particularly in 2°S–2°N. The SEC shear could be strengthened or weakened by the zonal wind, which results in higher or lower likelihood of MI and turbulence in the ELs. Moreover, tropical instability waves (TIWs) are energetic processes in the central-to-eastern equatorial Pacific, with the Yanai wave-initiated mode covering 4°S–4°N and the Rossby wave-initiated mode covering 5°S-8°N, both of which are associated with strong zonal velocity component49, providing additional vertical shear of zonal velocity22; in addition, the Yanai wave mode can also provide weaker stratification in particular phases, providing likelihood of MI from the stratification perspective50. Finally, the occurrence rate of EL mixing \({\gamma }_{{ELM}}\) in this region is prominent, with an average over 40%, confirming the feature of MI. Whereas, how well the schemes are suitable out of the 3°S–3°N band is hard to analyze, so that the estimation schemes are not adopted there.

Even with those uncertainties or lack of evidences, the results do show consistency with limited observations at 140°W, 0°N in the perspectives of daily, seasonal and interannual time scales, demonstrating that the present estimations are correct at least at the leading order:

  1. (1)

    At the daily time scale, the comparison between \({\varepsilon }_{{ELM}}\) and \({\varepsilon }_{{obs}}\) is shown in Fig. 2e, which shows that the most sample biases are within a factor of 1.5.

  2. (2)

    At the seasonal time scale, the comparison between estimated \({J}_{q}^{{EL}}\), its difference with \({J}_{q}^{s}\) (see below), and the resulted diabatic warming/cooling effect around 140°W, 0°N, with the observed parameters at 140°W, 0°N are shown in Supplementary Fig. 4, which show a significant controlling effect to the seasonal SST variation as revealed by ref. 31.

  3. (3)

    At the interannual time scale, the diabatic warming/cooling effect at the CP region (170°E–140°W) is shown in Fig. 4d, which resembles that at 140°W, 0°N during a different period (Fig. 4 of ref. 20).

Sampling error analysis on \({\varepsilon }_{{ELM}}\), \({k}_{{EL}}\), and \({J}_{q}^{{EL}}\)

One of the advantages of this work is to take full use of the long-term and wide coverage of Argo profiles, which maximizes the samples of subsurface turbulence data. In the analysis, we further accumulate the samples into large temporal and/or spatial bins and make mean of the properties, such as simultaneously over a long time (2000–2022) and a 5° × 6° box (for example, Fig. 1c), over a month and a large region, i.e., the CP (170°E–140°W, 3°S–3°N) or EP (130°–110°W, 3°S–3°N) (for example, Fig. 4a, e), and over several months of the growing El Niño /La Niña phase and CP/EP region (for example, Fig. 6c, d), where the samples are sufficient. The sampling errors for the sample-mean of the above properties are all estimated based on the 95% bootstrap CIs, and are shown in corresponding figures. We see that the mean values vary little along the equator (Fig. 1c), indicating small sampling errors; whereas, for the monthly means of earlier years (earlier than 2012), the samples are less, the mean values vary in a larger magnitude (Fig. 4a, e), but they do not change the conclusions of the present study.

Calculation of the heat flux retained in the ML from sea surface (\({J}_{q}^{s}\))

For each Argo profile, \({J}_{q}^{s}\) is calculated as \({J}_{q}^{s}={J}_{q}^{0}-{I}^{{h}_{{ELM}}}\), where \({J}_{q}^{0}\) is the net surface heat flux, and \({I}^{{h}_{{ELM}}}\) is the penetrative shortwave radiation that leaves the ML (at \({h}_{{ELM}}\)) further downward12,31, which is calculated as \({I}^{{h}_{{ELM}}}={C}_{{pen}}{I}^{0}\), where \({I}^{0}\) is the short wave solar radiation, and \({{C}_{{pen}}=0.62e}^{-{h}_{{ELM}}/1.5}+0.38{e}^{-{h}_{{ELM}}/20}\) is an empirical penetration coefficient model51. Both \({I}^{0}\) and \({J}_{q}^{0}\) are obtained from the daily ERA5 data, and are interpolated onto the positions and dates of the Argo profiles. The estimated \({J}_{q}^{s}\) also shows a perfect annual cycle, and jointly with \({J}_{q}^{{EL}}\), well displaying their controlling effect in the seasonal cycle of SST in the central equatorial Pacific Ocean (Supplementary Fig. 4; ref. 31).

Calculation of the contributing terms to CP and EP ENSO’s growth

The ML temperature equation is expressed as:

$$\frac{\partial {SST}}{\partial t}=-\underbrace{{{\vec{u}\cdot {\nabla }_{H}\left({SST}\right)}}_{{adv}}}+\overbrace{\underbrace{{\frac{1}{{h}_{{ELM}}\left({J}_{q}^{s}\right.}}_{\frac{\partial T}{\partial t}|{J}_{q}^{s}}} {-\underbrace{{\left.{J}_{q}^{{El}}\right)}}_{\frac{\partial T}{\partial t}|{J}_{q}^{{El}}}}}^{\frac{\partial T}{\partial t}|\Delta {J}_{q}}-\underbrace{w\frac{\partial T}{\partial z}+\nabla _{H}(k_{H}\nabla _{H}SST)}_{{Resi}}.$$
(20)

Here, the first term on the rhs is the horizontal advection term (adv), where \(\vec{u}\) is the mixed layer velocity vector, obtained from the OSCAR dataset, and SST is obtained from the GHRSST dataset (see above); the second term on the rhs is the diabatic warming/cooling term (\(\frac{\partial T}{\partial t}|\Delta {J}_{q}\)), which is constituted with both the \(\frac{\partial T}{\partial t}|{J}_{q}^{s}\) and \(\frac{\partial T}{\partial t}|{J}_{q}^{{El}}\) terms; because \({h}_{{ELM}}\) is almost equivalent to \({h}_{{MLB}}\) (Supplementary Fig. 2b), this term-induced change rate of temperature over the \({h}_{{ELM}}\) can be considered as it-induced change rate of SST. The third and fourth terms on the rhs represent the vertical advection and sub-grid horizontal diffusion terms, respectively, which are unable to estimate due to lack of data; however, since the overall SST change rate term on the lhs can be obtained from the GHRSST dataset, they can be calculated as the residual term (Resi).

Monthly and area-mean of those terms are calculated over the period 2000–2022 and over the CP (170°E–140°W) and EP (130°–110°W) regions; subsequently, monthly climatology and monthly anomalies are obtained for further analysis; particularly, they are assessed on the growing months (SON) of El Niño and La Niña phases, which are shown in Supplementary Fig. 7.