Subsurface ocean turbulent mixing enhances central Pacific ENSO

Domain: earth_science
Published: 2025-03-08
Editorial outcome: major_revision

Abstract
Sincethe21stcentury,theElNiñoandSouthernOscillation(ENSO)exhibitsmorepronouncedsignalsinthecentralPacific(CP)ratherthantheeasternPacific(EP),butthepredictionskillhaswaned,suggestinglimitedunderstandingofcrucialdynamicswithinthepredictionframework.Theoceanmixingaroundthemixedlayerbase,whichtransfersheatdownwardinadiabaticmanner,wasconsideredapotentialinfluencingfactor;yet,itseffecthasnotbeenadequatelyexaminedineitherCPorEPregionsduetoinsufficientdata.Here,weproposeanArgoprofiledata-basedmixingestimationmodel,whichyieldsabundantestimatesofsubsurfaceoceanmixingandturbulentheatflux.Consequently,wefindsignificantpositivefeedbackoftheoceanmixing-induceddiabaticwarming/coolingontheCPENSO,butnotontheEPENSO.Particularly,thediabaticeffectdominatesseasurfacetemperaturechangeintheCPregion,highlightingthenecessityfordiabaticCPENSOpositivefeedbackdynamicsinpredictionmodels.

Introduction
ENSOpredictionskillhasnotshownobviousimprovementoverthepastfewdecades,andinfact,therehasbeenadeclineattheturnofthe21stcentury1.ThedecreasedskillislargelyattributedtotheincreasedoccurrenceofCPElNiños2–4. Previous theories of ENSO primarily focused on the positive Bjerknes feedback5 and several negative feedback6–11. These theories highlight the wind-driven redistribution mechanismofthewarmwaters,eitherbetweentheeasternandwesternequatorialPacific(knownasthezonaladvectionmechanism),orbetween the equatorial and off-equatorial regions (known as the thermoclinemechanism).Thesemechanismsareallessentiallyadiabaticfromtheoceanicperspective.

However,insightsfrombothnumericalmodelsanddiagnosticsoftheparameterizedverticaldiffusion/entrainmentoftemperaturesuggestthatdiabaticprocessesaroundthemixedlayer(ML)base(MLB)mayplayanimportantroleinthedevelopmentofbothCPandEPENSO12–17. Direct turbulence observations demonstrate that entrainment and turbulent mixing just beneath the MLB frequently occur fromthe western18 to centraland eastern equatorialPacific19. Recently, using approximately 10 years of nearly continuous turbulenceobservationsbymooredsensorsat140°W,0°(Fig.1a),WarnerandMoum20foundpositivefeedbackofsubsurfaceturbulentmixingatthissitetotheENSO-relatedinterannualseasurfacetemperature(SST)changes(representedbytheNiño3.4region-basedOceanicNiñoIndex).Whilerepresentingagreatadvanceindata-baseddiagnosticofturbulentmixing’seffectonENSO,theemployedobservationsiteisunfortunately notlocated inthe core area ofeither CP ENSO or EPENSO (represented by Niño 4 and Niño 3 regions, respectively). ExtensivesamplingthatencompassesatleastthecentraltoeasternPacificequatorialbandarerequired,therefore,tovalidatethispositivefeedbackmechanism,andtodistinguishitseffectonEPfromthatonCP 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.

Results
Theentrainmentlayermixingontheequator
We firstfocus on identification of mixing events within a layer surroundingtheMLB.ThereasonwhywefocusonsuchalayerisbecauseitisjustbeneaththeMLBthatturbulentmixingdirectlytransferstheheatintheMLdownwardandhenceinfluencestheMLtemperature.Thislayerisreferredtoastheentrainmentlayer(EL)henceforth(seeMethodsfordetaileddefinitions).Themixingeventsareidentifiedasdensityoverturns21–23, because the latter is one of the key features duringshearinstability-associatedbillowsandturbulencesintheshearflow21.ThemixingeventsaredeterminedacrossthelowlatitudePacificOcean from the long-term (2000–2022), wide-coverage, high-resolution and quality-controlled Argo profile dataset; there, the Argoprofilesareconcentratedtotheequatorialband(Fig.1b),withmore than 40,000 Argo profiles available in the central-to-eastern equatorialPacificOcean(170°E–100°W,3°S–3°N),providingabundantdatasources.AteachArgoprofile,amixingeventisdeterminedonlywhen the density overturn patch covers three or more consecutive verticalgridpoints,andanELmixingeventisdeterminedonlywhenpartorfullofadensityoverturnislyingwithintheoutmostboundsofthe EL (see Methods). Combined with a velocity shear-based, wind-driventurbulentmixingscalingschemeasdescribedbelow,theArgoprofilesinthisregionleadtoabundantsubsurfaceturbulentheatfluxestimates;further,associatedwiththewidelyavailablesurfacenetheatfluxes(Fig.1a,c),thediabaticeffectontheMLtemperaturechangecanbeobtainedeventually.

Beforegivingthedetailsoftheturbulentmixingscalingscheme,wefirstshowtheoccurrenceoftheELmixingeventsandanalyzetheirformationmechanisms.ThedetecteddensityoverturnsrevealwidespreadandextensiveemergenceofELmixingeventsthroughouttheequatorialPacificOcean(Fig.2a),withaconcentrationofoccurrencerate (r ) in the central-to-easternequatorialPacific (170°E–110°W, 3°S–3°N),amaximumr of>50%centeredat160°W(Fig.2b).Themeanr overthisconcentrationregion(170°E–110°W,3°S–3°N)is41.2%;here,r ismeridionallyconsistentfrom3°Sto3°S(Fig.2c).Incontrast,r diminishesto<20%westof170°Eandeastof110°W.Between140°–170°E(thePacificwarmpool),r isalsomeridionallyconsistent; 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) experiencesmorefrequentdiapycnalcommunicationsbetweenthesurfaceandsubsurfacelayersthanboththefareasternPacific(eastof110°W)andthewesternPacific(WP,westof170°E)regions.

Firstly,windstressforcingoffersoneofthebasicinterpretationsforoccurrenceofELmixing.Strongerzonalwindstress(strongertradewind)leadstohigheroccurrenceofELmixingviatwomechanisms.Thefirstisitprovidesstrongerverticalshearandhencemoreshearinstabilityconditioninsubsurfacelayer,whichwillbedemonstratedinthefollowingsubsection,andthesecondisitleadstostrongerandquicker enhancement of entrainment turbulence around the MLB, whicharedemonstratedbylaboratoryexperiments24,directnumericalsimulation25, and observed by several years of turbulence measurements19.Ourresultsshowlarger andobviousresponseof r to wind stress at multiple time scales in the central-to-eastern equatorialPacific,because,inthisregion,boththemeanmagnitudeand variations of the trade wind are strong. In contrast, the WP (<160°E)ischaracterizedbynear-zerowindsandhencefewmixingevents(SupplementaryFig.1h).Inthesoutheastoftheequator,themeanzonalwindstressislargerthannortheastoftheequator,whichisassociatedwithhigherr there.Thewind’seffectonr isparticularlyobviousintimeseries.Inthehighr ELM region(170°E–110°W, 3°S–3°N),r exhibitsdistinctseasonalvariation and ENSO-related temporal scales. To emphasize the ENSO-focused interannual variations,13-monthrunningmeansofboththemonthlyr andτ arecalculated, which show a statistically significant correlation (with a peakcorrelationcoefficientR=–0.75,95%CI=[–0.82–0.64],andP-value<0.001over2010–2022).

Secondly, local background shear stability condition26 offers anotherbasicinterpretationforoccurrenceofELmixing.TheoccurrenceofELmixingreflectsfulfillmentoflocalshearinstabilitytherein.On average, the subsurface of the central equatorial Pacific (170°E–140°W) features as moderate background vertical shear of zonalvelocityfromtheSouthEquatorialCurrent(SEC)andtheEquatorialUndercurrent(EUC),andweakstratificationduetothePacificwarmpool,whilethenear-easternPacific(140°–110°W)isfeaturedasstrongshearandmoderatedstratificationduetotheshoalingofthethermocline (Supplementary Fig. 1a–f); as a result, this central-to-easternequatorialPacificregion(170°E–110°W)hasmorelikelihoodtoformastateofmarginalinstability(MI,i.e.,theRichardsonnumberRi≈0.25;SupplementaryFig.1g),andresultsinfrequentintermittentmixingevents27,28.ThecenteroflowRilocatesat160°W,thesameasthecenterofhighestr .Incontrast,inthefareasternPacific(eastof110°W),thoughtheshearofzonalvelocityisstrong,thethermoclineisverysharpandhencethestratificationisstrong,whileinthewesternPacific(westof170°E),thoughthestratificationismostlymoderate,theshearistooweak;asaresult,thelikelihoodofshearinstabilityinbothregionsisgreatlyreduced,soasthemixingevents.Therefore,thezonalpatternofsubsurfaceRi(SupplementaryFig.1g)mimicsthatofr (Fig.2b).

ItisnotedthattherolesofwindstressandshearinpromotingthegenerationofELmixingarenotindependent,particularlyunderMIconditionwhichisfrequentlyfulfilledinthecentral-to-easternequatorial Pacific region (170°E–110°W). On the one hand, wind stress controlsthestrengthofsubsurfacebackgroundshearandtheinitiationofentrainmentattheMLB,ontheotherhand,strongbackgroundshear or MI condition favors shear instability and the wind-driven entrainmenttodevelopintodeepcycleturbulence.Overall,theMIinthesubsurfacelayerprovidesthepreconditionfortheoccurrenceofturbulence,whilethewindstressinitiatestheturbulencethere.

Accordingly,wedevelopamethodtoquantifythemixingcoefficients,includingturbulentkineticenergydissipationrate,turbulentdiffusivityandthusturbulentheatfluxes,basedonbothwindstressandlocalshear.Thisenablesustorevealthespatial-temporallyvaryingELmixingacrossboththecentralandeasternPacificregions.WestartfromarecentworkbyWhittetal.29,whichproposedascalingmodel(hereafterWhittetal.’sscalingmodel)forthemaximumverticalturbulentbuoyancyfluxes(F )ofstrongmixinginthesubsurfacelayer,expressed as F (cid:1)0:2u2S, where S is the vertical shear of the horizontal velocity, and u2 represents the wind stress-induced friction velocity squared(u2=jτ jρ(cid:3)1),withρdenotingwaterdensity.Appro* x priately,thedetectedELmixingevents(overturns)inthepresentstudycanberegardedasthe‘strongmixing’neartheMLB,whichinducesmaximum turbulence buoyancy/heat fluxes therein. Furthermore, employingtheturbulentkineticenergyequation,andassumingthatthe turbulence is undergoing a nearly steady state (Ri≈0.25, i.e., S(cid:1)2N)andtheturbulentkineticenergyisonnearlyequilibrium30,weextendWhittetal.’sscalingmodelandproposeestimationmodelsforboththeturbulentkineticenergydissipationrate(ε )anddiapycnaldiffusivity (k ) of the EL mixing in the equator (within ±3°, 160°E–110°W),i.e.,

ε (cid:1)1:6u2N, ð1Þ
ELM *

and

k (cid:1)0:32u2 * =N, ð2Þ
ELM

whereN representsthemeanbuoyancefrequencyovertheEL.Surprisingly,theschemesshowgoodconsistencywithobservationaldataat140°W,0°N(Fig.2e;refertoMethodsfordetailed derivationsofEqs.1and2,aswellasuncertaintyandvalidationanalysisfortheseestimationschemes).Theadvantageoftheseschemesliesinthat,theirdeterminativefactorsu2andN,aswellasthestrongmixingcondition,canbedirectlyobtainedfromavailablewindstressandArgoprofiledatasets.

Conversely, for weakornon-ELmixing conditions whereWhittetal.’sscalingmodelisinapplicable,awidely-usedbackgrounddiffusivityk =10(cid:3)5m2s(cid:3)1isprescribed(Methods).

Finally,wepresenttheELmixingdiffusivityk foralltheArgoprofilesintheequatorialregionbycombingthetwo conditions:

k = { 0:32u2=N, strongELmixingexists
EL 10(cid:3)5(m2s(cid:3)1), noELmixingdetected

Theobtainedmeank EL (includingbothk ELM andk ELnomix )alongtheequatorialbandisshowninFig.2f.Itexhibitsanaveragemagnitudeof10−4 m2 s−1,yetisseveraltimeshigherin theverycentralequatorialPacificregion(180°–140°W)comparedtoeitherwestoreastofit.

Thedeterminationofk providesthebasisfordirectestimationofturbulentheatfluxinthe E E L L,whichiscalculatedas:

JEL=ρC k T , ð4Þ
q p EL z

whereρiswaterdensity,C p theheatcapacity,andT z thebulkmeanverticalgradientofsorted temperatureaveraged either over theELmixingpatchesforthestrongELmixingcategory,orovertheELforthenon-ELmixingcategory(Methods;SupplementaryFig.2a).Here,JEL representsthediapycnalturbulentheatfluxdownwardfromthemeandepthoftheEL,h ,whichcloselyapproximatesthedepthoftheMLB(h )basedonstatistics(SupplementaryFig.2b).

Meanwhile, at each Argo profile location, the net inward diapycnalheatfluxfromtheseasurfaceintotheML,Js,iscalculatedbysubtractingaportionofshortwaveradiationthatpenetratesfurtherdownward at h (IhELM) from the net surface heat flux (J0) (Methods).

As a result, long-term, wide-coverage data of JEL and Js are quantified(SupplementaryFig.3a,b),whichiscrucialforassessingthediabaticcontribution to temperature variations in the ML across varioustimescalesandregions.ThetemporalmeansofJE q LandJs q alongtheequatorialband(within±3°)showcomparablemagnitudes,rangingfrom~10Wm−2at170°Eto~150Wm−2at110°W,andprominentseasonalandinterannualvariations(alsoseeFig.1c).BothJEL andJs exhibitmagnitudesandseasonalvariationsquitelikethoseobservedat 140°W, 0°N (Supplementary Fig. 4a; ref. 31). The magnitudes of ENSO-focusedinterannualvariationsforbothJE LandJs reach50Wm−2(Fig.1c;SupplementaryFig.3c,d),indicatingthatJELcanexertagreatinfluence on ML temperature comparable to that of Js q at the ENSO scale.

Theentrainmentlayermixing’sresponsetoENSO
Horizontally,themonthlyanomaliesofJs (i.e.,Js0 )showaclearout-of-phase relationship with the ENSO-related interannual varying SSTanomalies (SSTAs), and demonstrate nearly zonal coherence acrossthe entire central-to-eastern Pacific basin (170°E–110°W) (SupplementaryFig.3c).Incontrast,themonthlyanomaliesofJEL (i.e.,JEL0 )displayanapproximatelyoppositecoherentpatternbetweentheCPregion(whichisdefinedspecificallyfor170°E–~140°Wand3°S–3°N)andtheEPregion(whichisdefinedspecificallyfor~130°–110°Wand3°S–3°N), although with their border migrating east- and westward year-by-year(SupplementaryFig.3d).

NotethatboththeCPandEPregionsdefinedinthepresentstudyaredifferentfromtheconventionalNiño4andNiño3regions(seeFig.2b,f),respectively.WhiletheNiño4andNiño3regionsrepresentthecoreregionsonSSTanomalyofCPENSOandEPENSO,respectively,thedefinedCPandEPregionsheresignifythezonaldiscrepancyofthediabaticeffectontheinterannualSSTchanges.Wewillshowthatthis pattern is associated with the zonal discrepancy of the ENSO-focusedinterannualvariationsofwindstressesinthefollowingsections.Zonally,theCPregioncoversmostoftheNiño4regionandalsothe western part of the Niño 3 region, so that it can represent the conventionalCPENSOregion,andcanrepresentthewesternpartofthe conventional EP ENSO region as well; in the meanwhile, the EP regionasdefinedinthepresentstudycanonlyrepresenttheeasternpartoftheconventionalEPENSOregion.

IntheCPregion,theeasterlytradewindsundergostrongintensificationduringthegrowingmonths(SeptembertoNovember,SON)oftheLaNiñaphase(Fig.3b,c;SupplementaryFig.5a),whichimpartgreater momentum to the upper ocean. Moreover, the enhanced easterlywindsstrengthenthewestwardSECintheupperlayer,andinducesenhancedeastwardEUCinthesubsurfacelayerbyresultinginwest-east pressure gradient32. Consequently, the squared shear (S2 u =ð∂ ∂ u z Þ2 ) of the zonal velocity (u) is amplified within the EL (and further below); meanwhile, the stratification (N2) remains almost unchanged inthe CP region, or even slightlydecreased particularlywestof160°W.Consequently,asdisplayedinSupplementaryFig.5c–f,thereducedshearsquared,S2 =S2(cid:3)4N2 (orthereciprocaloftheRichardsonnumber,Ri (cid:3)1= S2 u) u , Re b d ecom u esanomalouslypositive,whichprovidesmorefavorableshe N a 2 rinstabilityconditionsandhencemixingevents (Fig. 3e, f). The wind anomalies reverse during the growingmonthsofElNiñophase,whichultimatelysuppressshearinstabilityandturbulentmixingtherein(Fig.3a,d,c,f;SupplementaryFig.5a,b).

Theanalysesabovearerepresentedbythedifferencesinu,S2,N2 andS2 betweenwestwardandeastwardanomalouswindstressesat170°W,0°N(SupplementaryFig.5c–f),wheredataisavailablethankstotheTAOproject(seeMethods).Thisresultalsocanbeinferredfrom140°W,0°N20.

Consequently, more EL mixing is detected during the growingmonthsofLaNiñathanthatofElNiñointheCPregion(Fig.3f;SupplementaryFig.5b).Therefore,ther (Fig.4b),themeank (Fig.4c)andthemeanJEL(Fig.4d)thereinallshowstrongcorrelationswiththeENSO-relatedwindstressτ .[Inaddition,timeseriesalsoshowsthestrongcorrelationbetweenmagnitudesofτ andJELatseasonaltimescales(Fig.4a).]TheseresultsclearlyillustratetheprojectionofENSO-driveninterannualvariabilityacrosstheair-seainterfaceontothesubsurfaceturbulentmixing.

ComparedtotheCPregion,theEPregionexhibitsdifferencesinbothwindpatternsandshearinstabilityconditionsduringtheENSOcycle. Firstly, the interannual variation of zonal wind stress is weaktherein(Figs.3a–c,4e),withslighteasterly(westerly)anomaliesduringthegrowingmonths(SON)ofElNiño(LaNiña)phase(Fig.3a–c),whicharewestward(eastward)anomalouswindstresses,contrarytothatoftheCPregion.TheinfluenceofthisonELmixingisnotobvious,withtheinflectionpointofinterannualvariationsinr ELM andτ x occurringat~135°W.Thewind-drivenshearsquaredS2 intheEPregionexperiencesaslightincrease(decrease)duetotheweakwestward(eastward)windanomaliesduringthegrowingmonths(SON)ofElNiño(LaNiña),followingasimilarmechanismasintheCPregion.Itfavors(damps)theshearinstabilityfromtheperspectiveofshearcomponent.However,thestratification(N2)intheEPregionbecomesweaker(stronger)duetodeepening(shoaling)ofthethermocline,resultinginslightlyfavorable(unfavorable)shearinstabilityconditionsfromtheperspectiveofstratification.Asaresult,comparedtotheCPregion,thereisanobviousbutweakincrease(decrease)inELmixing(includingr ELM and the mean k ; Fig. 4f, g) in the EP region during the growingmonths of El Niño (La Niña), which reflects the primary drivingmechanismofwindstress(Fig.4e).Theinflectionpointthatseparatesthepositiveandnegativer interannualanomaliesalongtheequatorlocatesalsoaround135°W E ( L F M ig.3f),closetotheinflectionpointforthepositiveandnegativewindstress’sinterannualanomalies(Fig.3c).TheinflectionpointdeterminestheborderbetweenCPandEPregionsasdefinedinthepresentstudy.However,likelyduetothereducedverticalgradient,theinterannualvariationofthemeanJELisdiscrete,thusitshowsalesssignificantrelationshipwithτ inthisregion(Fig.4h).

Theentrainmentlayermixing’sfeedbackonENSO
Now, we examine the diabatic effect on ENSO. A net downward Js warms the ML water, while anetdownward JEL cools it.Jointly,thedifferencebetweenJs andJEL,i.e.,

ΔJ =Js (cid:3)JEL, ð5Þ
q q q

representstheoverallretainedheatviadiapycnalheatfluxesintheML,withthediabaticwarming/coolingratebeingestimatedas20,31:

(cid:2) (cid:2) T t (cid:2) (cid:2) (cid:2) Jq = ∂ ∂ T t (cid:2) (cid:2) (cid:2) Jq = ρ 0 C Δ p J h q ELM : ð6Þ

Similarly, a positive (negative) monthly anomalous Js0 results in anomalousdiabaticwarming(cooling)oftheMLwater,whileapositive(negative)monthlyanomalousJEL willleadtoanomalousdiabaticcooling(warming).Meanwhile,theirdifferences,denotedasΔJ 0andJEL0 ,respectively,representtheoverallanomalousdiabaticwarming/coolingeffectsontheMLwater.

InadditiontotheoutstandingannualcycleofΔJ ,whichgovernstheseasonalvariationoftheSSTinboththeCPandEPregions(Fig.5a;SupplementaryFig.4a–c;ref.31),theinterannualvariationinΔJ 0also standsout(Fig.5b;SupplementaryFig.3c,d),whichshowpotentialrelationtotheENSOcycles(Fig.5c)

Wefirstexplorethequalitativediabaticwarming/coolingeffectonthegrowingofElNiñoandLaNiñainboththeCPandEPregions.Onaverage,duringthegrowingmonths(SON)ofElNiñophases(Fig.6a),both JEL 0 and Js0 are negative in the CP, however, Js0 is smaller inmagnitude compared to JEL, which thus jointly lead to anomalousdiabatic warming (ΔJ 0>0) of the ML water. Similarly, during thegrowingmonths(SON q )ofLaNiñaphases(Fig.6b),bothJEL 0 andJs0 arepositive, but the former outweighs the latter, jointly leading toanomalous diabatic cooling (ΔJ 0<0) of the ML water. The resultsdemonstrate the positive role of the diabatic warming/coolingmechanisminpromotingthedevelopmentofbothElNiñoandLaNiñaintheCPregion.ThispositiveroleisprimarilyinducedbyJEL.

Inthemeanwhile,intheEPregion,Js0 issignificantlyneg q ativewithalargemagnitudeduringthegrowingmonths(SON)ofElNiñophase,leading to anomalous cooling of the ML water therein, while JEL isweakandpositive(significantonlyeastof120°W),alsocontributingtoanomalous warming. Together, they induce significant anomalousdiabatic cooling (ΔJ 0<0) (Fig. 6a). In contrast, during the growingmonthsofLaNiñap q hase,thestrongpositiveJs0 andweakJEL0 leadtoanomalouswarming(ΔJ 0>0)therein(Fig.6b).ItisobviousthatthisoveralldiabaticeffectplaystorepressthegrowthofboththeElNiñoandLaNiñainthisEPregion,andiscontrolledbytheJs0 component.

Nowlet’squantifytheircontributionstotheSSTchange.IntheCPregion(Fig.6c),duringthegrowingmonthsofElNiño,themeanSSTAincreasesatarateof~+0.2°C/month(alsoseeSupplementaryFig.6g),while the diabatic warming mechanism contributes to ~+0.35 °C/month(alsoseeSupplementaryFig.6e),obviouslyoutweighingotherfactors and dominating the warming trend. Similarly, during thegrowing months of La Niña, the mean SSTA changes at a rate of~–0.14°C/month), whereas the diabatic cooling mechanism contributesto~–0.29°C/month,onceagaindominatingthecoolingtrend.Notably, these contributions primarily stem from the larger ENSO-relatedJEL0 ,ratherthanJs0 (Fig.6c;SupplementaryFig.6a,c,e),whichisassociatedwithratesof+0.48and–0.15(–0.50and+0.22)°C/month,respectively,forthedevelopmentofElNiño(LaNiña).Moreover,thepositivecontributionfromthediabaticwarming/coolingmechanismexceedsthepositivecontributionfromthehorizontaladvectionintheCP region, both of which should be balanced by the negative contributionfromtheverticaladvectionandhorizontaldiffusionterms(SupplementaryFig.7a,b).

Incontrast,thediabaticwarming/coolingmechanismexhibitsanegativeeffectonthedevelopmentofeithertheElNiñoorLaNiñaintheEPregion(Fig.4b,c).DuringthegrowingmonthsofElNiño,theMLwater warms at a rate of ~0.21°C/month, while Js0 and JEL0 lead toanomalouscoolingatarateof~–1.0°C/monthandwarmingatarateof~+0.24°C/month, respectively. Obviously, while JEL0 promotes thegrowth of El Niño, Js0 plays a dominant damping role; in sum, thediabatic mechanism represses the growth of the El Niño therein.Similarly,duringthegrowingmonthsofLaNiña,theMLwatercoolsat a rate of ~–0.21°C/month, while JEL0 leads to anomalous cooling of~–0.18°C/month, Js0 leads to diabatic warming of ~+1.5°C/month,dominantthedampingoftheLaNiñahere(alsoseeSupplementaryFig.6f,h).

ThisdiscrepancyrelativetotheCPregionarisesduetothatΔJ0isdominatedbyJs0 inthisregion(Fig.6a,b,d;SupplementaryFig.6b,d,f), which varies inter-annually out-of-phase with SSTA, while JEL0 variesweakly,discretealongwithwindstress(Fig.4h),andonlysignificant east of 120°W, hindering its diabatic effect on the ENSOdevelopmentintheEPregion.Inotherwords,thediabaticwarming/coolingmechanismactstosuppressthedevelopmentofENSOintheEP region (Supplementary Fig. 6f); this negative contribution isreversed by the positive contribution fromboth the horizontal andverticaladvectioninthisregion(SupplementaryFig.7c,d).

Insummary,thediabaticwarming/coolingmechanismplaysdifferent roles in the two regions. In the CP region, it serves asaprominentpositivemechanismtothedevelopmentofbothElNiño orLaNiña.Inturn,thispositivemechanismisintrinsicallytiedtotheELmixing,whichisfurtherprimarilygovernedbytheENSO-associatedwind stress variations. Whereas, in the EP region, it represses thedevelopmentofENSOprimarilythroughtheJs0 component.

Methods
TheArgoprofiledata
TheArgo(ArrayforReal-timeGeostrophicOceanography)profiledata 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 2000mintheocean,withanaccuracyof±0.002°Cfortemperature,± 0.005 psu for salinity, and ± 2 decibar for pressure35. Only those quality-controlledprofileswithaverticalresolutionof2morlessareused. Totally >160,000 Argo qualified profiles are available in the tropicalPacificOcean(Fig.1b).

TheERA5reanalysisdata
Thezonalwindstress(τ ),netsurfaceheatflux(J0)andshort-wavesolarradiation(I0)componentsofthisdatasetareusedfortheestimationofJELandJs(seebelow).TheERA5dataisthefifthgenerationofECMWFreanalysis36,coveringtheperiodof1940topresent.ERA5isofhigh spatial (0.25°) and temporal (hourly) resolutions. An obviousimprovement is observed in the ERA5 over previous versions. Thisdataset is available at https://doi.org/10.24381/cds.adbb2d47. Daily meansareusedforanalysis,andthevariablesareprocessedontothepositionsanddatesoftheArgoprofiles.

Thelong-termTAOdata
TAO is the abbreviation for the Tropical Atmosphere and Ocean project37.Thisdatasetprovidesoceanvelocity,temperature,salinityprofiledataandwindstressinthecentraltoeasterntropicalPacific,available at http://www.pmel.noaa.gov/tao/data_deliv. Daily ocean data at 170°W, 0°N and wind stressdata in the vicinity are used in determination of conditions of anomalous westward and eastward windsnear170°W,0°N.Theanomalouswestwardandeastwardzonalwindconditions are obtained during 2008–2020 and over 180°–160°Wand5°S–5°N(a20°×10°regioncenteredat170°W,0°N).

TheIAPdata
Thisdatasetconsistsofmonthlytemperatureandsalinitydatafrom1960topresent,withahorizontalresolutionof1°×1°and41verticall evelsfromthesurfaceto2000m38.ItisusedforcalculatingtheN2anomaliesinFig.7.Thedataisfreelydownloadedfrom:http://www.ocean.iap.ac.cn/.

TheGHRSSTdata
TheLevel4ofGroupforHighResolutionSST(GHRSST;https://www.ghrsst.org/)39isadoptedforcalculatingtheSSTchangerateduringElNiñoand LaNiña phasesinboth theCP and EPregions (Fig.6c,d;SupplementaryFig.6g,h).TheL4griddedproductsaregeneratedbycombiningcomplementary satellite and in situ observations withinOptimalInterpolationsystems.ThisdatasetprovidesdailySSTat0.25°horizontalresolution,andisdownloadedfrom:ftp://ftp-oceans.ncei.noaa.gov/pub/data.nodc/ghrsst/L4/GLOB/JPL/MUR/.

TheTIWE
TheTIWEdata
ThisdatareferstotheturbulencemeasurementsduringtheTropicalInstabilityWaveExperiment(TIWE)duringthefallof1991at140°W,0°N40.Nearly6000castsofmicrostructuretemperature,conductivity,and shear measurements in the upper 200m were made. Simultaneoushorizontalvelocitywasmeasuredbytheship-mountedAcousticDopplerCurrentProfilers(ADCPs).ThedatawasprovidedbyDr.Lienviapersonalcommunication. The data were averaged daily into 36 profiles with the vertical resolution of 1m from November 4 tillDecember12.

TheJohnsonetal.(2002)data
This dataset refers to the potential temperature, salinity, potentialdensity,andvelocitydatathatareobtainedintheequatorialband41,whichisdownloadedfromhttps://floats.pmel.noaa.gov/sites/default/files/atoms/files/meanfit_m.cdf_.zipviaDr.Johnson’s website https://floats.pmel.noaa.gov/gregory-c-johnson-home-page.

TheOSCARdata
ThisdatasetreferstotheOceanSurfaceCurrentAnalysesReal-time(OSCAR), which contain ocean surface mixed layer velocities calculatedfromsatellite-sensedseasurfaceheightgradients,oceanvectorwinds, and SST fields using geostrophy, Ekman, and thermal winddynamics42.Thev2isemployed,whichhasahorizontalresolutionof1/3°andtemporalresolutionof5daysfortheperiodof1993–present.Itis downloaded from: http://apdrc.soest.hawaii.edu/erddap/griddap/hawaii_soest_bb96_63ae_8950.html.

Definitionsofthemixedlayer(ML)andentrainmentlayer(EL)ofArgoprofiledata
The ocean surface ML is a physically well-defined layer and can be determinedbymanypreviouslyproposedcriteria.Inthisstudy,theMLisdeterminedsuchthatthedensityatitsbaseis0.01kg/m3greaterthanthenear-surface(at10mdepth)43.TheELisalayerthatcoversboththelowerflankoftheMLandtheupperflankofthepycnocline,so thatinternalwavesandshearinstabilitycanexist,whichthusinduceentrainment, turbulence and finally diapycnal heat transport across theMLbase(MLB).However,nocommonlyapplicablecriteriawereavailableyetforitsdefinition.Inthisstudy,forthesakeofobjectivity,wedefinealayerof20mthickastheEL,withitstopbeing5maboveanditsbottombeing15mbelowtheMLB(SupplementaryFig.2a).Thisdefinitionisachievedalsobyconsideringthefollowingreasons.Firstly,thepurposeofdetermininganEListodetectturbulentmixingthatmayinducecross-MLBturbulence,sothattheverticalboundariesofELshouldnotbefarawayfromtheMLB.Secondly,observationsshowthat, in either the EP29 or western equatorial Pacific18,23, maximum mixing coefficient or turbulent buoyancy/heat fluxes often occur within30mbelowtheMLB;whileofftheequator(at3°N),largeeddysimulation(LES)withasrealisticaspossibleoceanicandatmosphericconditions shows that maximum heat flux often occur within 10mabovetheMLB29.Giventhoseobservationalevidences,andconsideringotherwisethatdistantmixingeventsmaynotimpacttheML,wechoosethehalfdistances(5mand15mrespectively)totheMLBastheoutmostupperandlowerboundsoftheEL,whichisschematizedinSupplementaryFig.2a.

DeterminationofmixingeventsintheELandthemeandepthofELmixing(h ELM)
Inthepresentstudy,weadoptthedensityinversion(ordensityoverturn)detectionmethod21tothelong-term(2000–2022),widelycovered,high-resolutionquality-controlledArgoprofiles(seeabove),todetectmixingeventsintheEL.AnELmixingeventisdeterminedonlywhen the overturnpatchcovers three or moreconsecutive verticalgridpoints,andpartorfullofitislyingwithintheoutmostboundsofEL(SupplementaryFig.2a).SincetheverticalresolutionoftheArgoprofiledataisabout1–2m,theminimumdensityinversionscalethatcan be resolved is 4m; as such, only the large-scale overturns aredetected,whichrepresent‘strongmixing’therein.Smallscalemixingeventspossiblyexistbutcan’tbedetected.Therefore,itisreasonabletoassumethatthedepthswithdetectedmixingeventsisassociatedwith T la h r e ge m st e t a u n rb d u e l p e t n h t o fl f ux E e L s m . ixingforeachArgoprofile,denotedas h ,shouldalsobedetermined.IfmultipledensityoverturnswereELM detected within the EL, the mean depth of the largest overturn isdefinedastheh ;ifallthedensityoverturnshavethesamesizes,theELM meandepthofallofthemisdefinedastheh .Whereas,whennoELM densityoverturnwasdetectedwithintheoutmostELbounds,theMLBisdefinedastheh (SupplementaryFig.2a).Themeanh alongELM ELM theequatorforthetwocategoriesareshowninSupplementaryFig.2b.

TheproposedscalingschemesforestimationoftheELmixingcoefficientsεandkbasedonWhittetal.’sscalingscheme
Intense EL turbulence is frequently generated by surface wind andbuoyancy forcing under sufficient shear conditions therein. Whittetal.29proposedascalingschemetorepresentdailymaximumturbulentbuoyancyandheatfluxesinthesubsurfacelayerintermsofbulkverticalshearandwindstressforstrongmixing(whosemaximumbuoyancy flux roughly is larger than 10−7.5m2s−3). Theoretically, theturbulent kinetic energy (K) in a shear layer just beneath the MLBevolves due to shear production (F S) and dissipation (ðεÞ plusm buoyancyflux(F )intheshearlayer,whiletheshear(S)isdrivenbytheb verticaldivergenceofmomentumfluxes,whichiswrittenas26,30:

∂K
=F S(cid:3)ε(cid:3)F , ð7Þ
∂t m b

∂S = 1 ½u2(cid:3)F ðhÞ(cid:4), ð8Þ
∂t Hh * m

wherehisthedepthofMLB,Histhethicknessforbulkverticalshearcalculation,Sisthebulkverticalshear,F istheverticalmomentumm flux,andu2=jτjisthefrictionvelocitysquared,withτthewindstress* ρ andρthewaterdensity.Basedontheories,observationalevidencesand large eddy simulation (LES) results, Whitt et al. found that thefollowingrelationship,

F (cid:1)0:2u2S, ð9Þ
b *

isaplausiblescalingmodelwhichcanexplainmostofthedailyvarianceofF atthedepthofmaximumF .Thismodeldirectlyconnects b b thesubsurfaceturbulentbuoyancyfluxandbothsurfacewindforcingandlocalvelocityshearcondition.

However,thescalingmodel(Eq.9)isdifficulttoapplytoArgoprofiledatabecausethelatterdoesnotcontaininformationofvelocityshear.Therefore,weconductfurtherderivationstomakethismodelapplicabletotheabundantArgoprofiles.Theory(Eqs.7and8)suggeststhatwhentheshear(S)andturbulentkineticenergy(K)areinanearly steady state30, which is a good assumption at the daily timescale,thenu2(cid:1)F ðhÞandε(cid:1)F S(cid:3)F ;additionally,adoptingEq.9,* m m b weobtain:

ε(cid:1)0:8u2S: ð10Þ
*

Assuch,thisscalingmodeldirectlyrelatestheturbulentkineticenergydissipation rate (ε) in the subsurface layer to the sea surface windstress(τ)andthesubsurfaceverticalvelocityshear(S).Thisscalingmodelisconsistentwithatheoreticalderivation30,andissimilartoa previousscalingmodelthatisderivedforathickerlayerandhencehasasmallercoefficient44.

Inaddition,notethatwhentheturbulentkineticenergy(K)isina steadystate,thegradientRichardsonnumber(Ri)overtheshearlayer shouldbearoundthecriticalvalue,i.e.,

Ri= N S2 2 (cid:1)0:25, ð11Þ

otherwise, K will either increase (when Ri<0:25) or decrease(Ri>0:25)45.AstateofRi(cid:1)0:25 iscalledmarginalinstability(MI46),whichisparticularlytrueintheupperlayerofthecentralequatorialPacific27,28. Here, N= (cid:3) g ∂ρ is the buoyancy frequency and denotes ρ 0 g ∂ r z stratification,withg the ρ 0 g ∂ r z avitationalcoefficient,ρ =1025kgm−3the0 referencedensity,andρthepotentialdensity.Inastratifiedstate,abackground vertical velocity shear plays as a necessary conditionforMI27.Equation11thusleadstoS(cid:1)2N,andEq.10isfinallytransformedto

ε(cid:1)1:6u2N: ð12Þ
*

Theestimationscheme(Eq.12)directlyrelatesεinthesubsurfacelayertothewindstressτandsubsurfacebuoyancyfrequencyN,whichcanbemoreeasilyobtainedfromcurrentavailabledatasets.

Furthermore,thediffusivity(k)caneasilybecalculatedas:

k=Γε=N2, ð13Þ

withawidelyusedconstantmixingefficiency,Γ=0.247.UsingschemeEq.12,wefurtherrewriteEq.13as:

k=0:32u2=N: ð14Þ
*

Thisscheme(Eq.14)buildstherelationshipbetweenthesubsurfacediffusivitykandbothτandN.

Notethat,accordingtoWhittetal.’sscalingmodel,theobtainedε andkshouldberepresentativeforthepositionofmaximumF ,whileb NandSarebulkaveragedoveracertainthicknessinthevicinityofthemaximum F . Therefore, both the depth of maximum F and theb b thickness for the calculation of the bulk mean need to be furtherdetermined.

Firstly, for the case that Argo profiles are with detected ELmixingevents,themixingeventsareassumedasstrongmixingandleads to maximum F therein, because otherwise mixing eventsb won’t be detected, let alone turbulent fluxes. In this case (i.e.,Category1,whenstrongELmixingexists,asshowninSupplementary Fig. 2a), the mean depth of EL mixing, i.e. h , is naturallyELM considered as the depth of the maximum F , and also of the cal-b culated ε and k, again because this position is where the largestoverturnsaredetected.ThethicknessofELmixing(H )isdefinedELM as a 10-m thick layer that centered at h ; but, if necessary, isELM further constrained by the outmost bounds of the EL as definedabove (Supplementary Fig. 2a). In this category, the diffusivity iscalculatedasthesameasEq.14,i.e.,

k =0:32u2=N, ð15Þ
ELM *

whereNiscalculatedfromthesorteddensityoverH .ELM ELM Accordingly,thedissipationratecanbecalculatedfollowingthesameasEq.12,i.e.,

ε =1:6u2N, ð16Þ
ELM *

andalsocanbecalculatedbyEq.10whenthesheardataisavailable,i.e.,

ε =0:8u2S: ð17Þ
ELM *

Thecalculatedε withEqs.16and17showsconsistencyandmatchesELM theobservationsat140°W,0°Nverywell(Fig.1d).

Whereas,itshouldbenotedthat,inmanytimes,thewindmight be weak, and/or the shear stability condition is not fulfilled forentrainmentandturbulentmixingtohappenneartheMLB,oronlysmallturbulenceoccurs.Insuchcase,mixingeventscan’tbedetected in the EL on the Argo profiles (Case 2, Supplementary Fig. 2a), thediffusivityisthenisassumedtobesmallandprescribedawidelyusedbackgroundvalue,i.e.,

k =10 (cid:3)5m2s(cid:3)1: ð18Þ
ELnomix

Here,anELisstilldefined(a10mthicklayercoveringtheMLB)forthecalculationofturbulentheatfluxes.

Finally, the diffusivity of the EL for each Argo profile can bedescribedas:

( 0:32u2=N, strongELmixingexists
k = * ð19Þ
EL 10 (cid:3)5(m2s(cid:3)1), noELmixingdetected

Inthepresentstudy,u2 iscalculatedfromthedailyERA5zonalwind* stress.Statisticsshowthat,subjecttomonthsandspaces,theoccurrenceofELmixingγ rangesbetween20%~60%,andthe90per-ELM centile of k ranges between 10−2.8~10−2m2s−1 with the median ofELM 10−2.5m2s−1(SupplementaryFig.8a),whichismorethan2orderslargerthank .Takingbothk andk intoaccount(i.e.,k ),theELnomix ELM ELnomix EL medianof them reduces to 10−4.0m2s−1. ε of the EL mixing casesELM shows a mean magnitude of 10 (cid:3)8 Wkg−1 (ranging from 10−8.5 to10−7.2Wkg−1), and is larger in the eastern than central and western equatorial Pacific. Those parameters show clear seasonal and interannualvariations (Supplementary Fig. 8b). Allthose values and featuresfallintotherangeoflimitedobservations.

We note that, the schemes Eqs. 17 and 19 are only applied to3°S–3°Nand170°E–110°W,wherethebackgroundverticalshearnear theMLBcommonlyexistsandhenceWhittetal.’sscheme(Eq.9)29andthe assumption of MI are valid to a large extent (SupplementaryFig.1g).

Estimationoftheturbulentheatflux(JE q L)atthemeandepthofELmixing
For each Argo profile, either with or without detected EL mixingevents,JELiscalculatedasJEL=ρC k T ,whereC istheheatcapacity,q q p EL z p T isthebulkmeantemperaturegradientaveragedoverH forthez ELM k case,andalsoovera10-mthicklayeraroundtheh (i.e.,theELM h ) for the k case (Supplementary Fig. 2b). The latter onMLB ELnomix averageistwoorderssmallerthantheformer.ThecalculatedJEL isq prescribedath .TheestimatedJELshowsaquitereasonableannualELM ELM cycle(SupplementaryFig.4a),indicatingthatthescalingschemesfork andJELarevalidated.EL q

Uncertaintiesandverificationoftheaboveε ELM andJE q L estimations
Weareawarethatseveraluncertaintiesexistintheestimationmodelsofbothε andJEL:
(1) Accordingtoseveralpreviousstudiesascitedabove,toavoidextreme fluctuations, the mixing coefficients are usually estimatedatdailytimescalesasdonehere,however,duetolackofdata,thesnapshotNandT obtainedfromArgoprofileareusedz as the daily proxies. The influence of this assumption isunknown, but it is expected that the adopted temporal andspatialaveragemaysmoothedoutthebias.
(2) Theseasurfaceforcingvariables(windstressandsurfaceheatfluxcomponents)areadoptedfromtheERA5dataset,whichhasthelongestandfinestcoverage,showsgoodperformanceintheequatorial region48, and represents one of the best availablereanalysis data. However, any remaining uncertainty of thisdatasetinthestudyregionmightprojecttotheestimationoftheproposedestimationmodels.
(3) TheextensionofWhittetal.’sscalingmodelfrom140°W,0°Nand140°W,3°Ntothebroaderequatorialbandisnotdirectevidence-based,butby reasonable deduction with indirectinformation. The Whitt et al.’s scaling model is based onobservationsat140°W,0°Nandlargeeddysimulation(LES)resultsatboth140°W,0°Nand140°W,3°N.Thesuccessofthescalingmodelreliesoncommonexistenceofverticalshearin theupper layer and marginal instability (MI) condition inthe subsurfacelayer.Thisconditionislikelytobefulfilledalonethenarrowequatorialband(3°S–3°N)andinthecentral-to-eastern region (170°E–110°W) due to common existence ofbothverticalvelocityshearandnot-extremely-strongstratification therein, but no direct observational or simulationevidence available to our knowledge. However, a lot ofinformation indicate that MI could be occasionally occur.Basedonlinearstabilityanalysis,Zhangetal.28demonstratedthatMIoccupiesalargeportionoftimeintheupperlayerinthiszonalrangeattheequator.Intheoff-equatorialregion,thewestwardflowingSEC,whichcoversameridionalrangeofabout5°Sto5°N,providesthebackgroundverticalshearin the upper ~50m41. The SEC is the primary current thatprovidestheshearattheMLBandELdepths,andisenhancedbytheEUCparticularlyin2°S–2°N.TheSECshearcouldbestrengthenedorweakenedbythezonalwind,whichresultsinhigher or lower likelihood of MI and turbulence in the ELs.Moreover, tropical instability waves (TIWs) are energeticprocessesinthecentral-to-easternequatorialPacific,withtheYanaiwave-initiatedmodecovering4°S–4°NandtheRossbywave-initiated mode covering 5°S-8°N, both of which areassociatedwithstrongzonalvelocitycomponent49,providingY ad an d a it i io w n a a v l e ve m rt o ic d a e l c sh an ea a r ls o o f z p o r n o a v l id v e el w oc e i a t k y e 22 r ;i s n tr a a d ti d fi i c t a io ti n o , n th in e st a r r a t t i i c fi u c l a a t r io p n h p as e e r s s , pe p c r t o iv v e i 5 d 0 i . n F g in l a i l k l e y l , ih t o he od oc o c f ur M re I nc f e ro r m ate th o e f adv ∂ ∂ T t jJs q ∂ ∂ T t jJE q l Resi

Evenwiththoseuncertaintiesorlackofevidences,theresultsdo showconsistencywithlimitedobservationsat140°W,0°Nintheperspectivesofdaily,seasonalandinterannualtimescales,demonstratingthatthepresentestimationsarecorrectatleastattheleadingorder:
(1) Atthedailytimescale,thecomparisonbetweenε andε isELM obs showninFig.2e,whichshowsthatthemostsamplebiasesarewithinafactorof1.5.
(2) Attheseasonaltimescale,thecomparisonbetweenestimatedJEL,itsdifferencewithJs (seebelow),andtheresulteddiabaticq q warming/coolingeffectaround140°W,0°N,withtheobservedparametersat140°W,0°NareshowninSupplementaryFig.4,whichshowasignificantcontrollingeffecttotheseasonalSSTvariationasrevealedbyref.31.
(3) At the interannual time scale, the diabatic warming/coolingeffectattheCPregion(170°E–140°W)isshowninFig.4d,whichresemblesthatat140°W,0°Nduringadifferentperiod(Fig.4ofref.20).

Samplingerroranalysisonε ELM,k EL,andJE q L
Oneoftheadvantagesofthisworkistotakefulluseofthelong-term and wide coverage of Argo profiles, which maximizes thesamplesofsubsurfaceturbulencedata.Intheanalysis,wefurtheraccumulatethesamplesintolargetemporaland/orspatialbinsandmakemeanoftheproperties,suchassimultaneouslyoveralongtime (2000–2022) and a 5°×6° box (for example, Fig. 1c),over a monthandalargeregion,i.e.,theCP(170°E–140°W,3°S–3°N)orEP(130°–110°W, 3°S–3°N) (for example, Fig. 4a, e), and over severalmonthsofthegrowingElNiño/LaNiñaphaseandCP/EPregion(forexample,Fig.6c,d),wherethesamplesaresufficient.Thesamplingerrorsforthesample-meanoftheabovepropertiesareallestimatedbasedonthe95% bootstrapCIs,and areshownincorrespondingfigures.Weseethatthemeanvalues varylittlealongtheequator(Fig.1c),indicatingsmallsamplingerrors;whereas,forthemonthlymeansofearlieryears(earlierthan2012),thesamplesareless,themeanvaluesvaryinalargermagnitude(Fig.4a,e),buttheydonotchangetheconclusionsofthepresentstudy.

CalculationoftheheatfluxretainedintheMLfromseasurface(Js )
ForeachArgoprofile,Js iscalculatedasJs=J0(cid:3)IhELM,whereJ0istheq q q q netsurfaceheatflux,andIhELM isthepenetrativeshortwaveradiationthatleavestheML(ath )furtherdownward12,31,whichiscalculatedELM as IhELM=C I0, where I0 is the short wave solar radiation, andpen C =0:62e (cid:3)hELM =1:5+0:38e(cid:3)hELM =20 isanempiricalpenetrationcoeffipen cientmodel51.BothI0andJ0areobtainedfromthedailyERA5data,andq areinterpolatedontothepositionsanddatesoftheArgoprofiles.TheestimatedJsalsoshowsaperfectannualcycle,andjointlywithJEL,wellq q displayingtheircontrollingeffectintheseasonalcycleofSSTinthecentralequatorialPacificOcean(SupplementaryFig.4;ref.31).

CalculationofthecontributingtermstoCPandEPENSO’sgrowth
TheMLtemperatureequationisexpressedas:

∂T = (cid:3)~u (cid:5) ∇H T + 1 h EL J s q (cid:3) J E q l z (cid:3) w ∂ T z + ∇ H (cid:5) ð k H ∇ H S S T Þ: ð20Þ

Here,thefirsttermontherhsisthehorizontaladvectionterm(adv),where~uisthemixedlayervelocityvector,obtainedfromtheOSCARdataset,andSSTisobtainedfromtheGHRSSTdataset(seeabove);thesecondtermontherhsisthediabaticwarming/coolingterm(∂TjΔJ ),∂t q whichisconstitutedwithboththe∂TjJs and∂TjJElterms;becauseh∂t q ∂t q ELM is almost equivalent to h (Supplementary Fig. 2b), this term-MLB inducedchangerateoftemperatureovertheh canbeconsideredELM asit-inducedchangerateofSST.Thethirdandfourthtermsontherhsrepresent the vertical advection and sub-grid horizontal diffusionterms,respectively,whichareunabletoestimateduetolackofdata;however, since the overall SST change rate term on the lhs can beobtained from the GHRSST dataset, they can be calculated as theresidualterm(Resi).

Monthly and area-mean of those terms arecalculated over theperiod 2000–2022 and over the CP (170°E–140°W) and EP(130°–110°W) regions; subsequently, monthly climatology andmonthlyanomaliesareobtainedforfurtheranalysis;particularly,theyare assessed on the growing months (SON) of El Niño and La Niña phases,whichareshowninSupplementaryFig.7.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript “Subsurface Ocean Turbulent Mixing Enhances Central Pacific ENSO as Positive Feedback” by Liu et al\ndemonstrates the impact of vertical turbulent mixing on sea surface temperature variability in the equatorial Pacific. The\nauthors use density overturns in quality controlled ARGO profiles to infer mixing events around the entertainment layer. The\nauthors also assert that there is a more pronounced effect of the mixing feedback on Central Pacific ENSO events (CP)\nversus Eastern Pacific (EP) events. The method is innovative and a good use of existing data from Argo to estimate a\nquantity that is elusive. There are some questions I believe should be addressed before the manuscript should be fully\nconsidered for publication.\nMajor comments\nThe authors rely on the state of marginal stability as given by Richardson number around 0.25 to substitute stratification for\nshear in their relationships they employ to estimate Jq. The authors rightly state in the extended section that the assumption\nof marginal instability only holds in the presence of strong background currents inducing vertical shear. In this manuscript\nthis is assumed to be true within 3S-3N around the equator. However, the observed undercurrent which gives rise to the\nstrong current shear is shown in for example Johnson et al (2002) or Roger (Karnasukas et al 2020) to span at most 2S-2N.\nA rigorous analysis of which profiles, at which location are analyzed needs to be presented, and reasonable arguments\nneed to be made why the marginal stability assumption holds, even away from the strongest current system. Whitt et al (as\nreferenced in the manuscript) examine 3N partially in relation to strong mixing induced by tropical instability waves, no\nmention of TIW kinetic energy is made here. In the introduction of the text the authors reference the number of Argo profiles\n5S-5N. These discrepancies need to be reconciled. Why are we shown the number of Argo profiles in the 3-5N/S bands, if\nthey are (rightly) not analyzed? It is also not very clear in the body of the text how the density overturns relate to the shear\ndriven mixing. This is explained in the extended data, but should be mentioned concisely in the main text as well. Related: In\nFigure 3 we are looking at the response of mixing indicators to wind stress, which is again not entirely clear to me why. It\nneeds to be better explained in the text how this relates to shear driven turbulence.\nThe discussion about the Central Pacific (CP) and Eastern Pacific (EP) is at times not clear. In the abstract it is said that the\nCP exhibits a stronger signal, what is meant? Historically EP ENSOs are stronger than CP ENSOs. Do the authors mean\nthat CP is more frequent in the past two decades? Further, the authors speak about CP and EP but do not show coincidence\nof mixing during EP or CP events, just the regions. Side note / related minor comment: in some figures it says that CP and\nEP regions are marked, but I cannot see the marking, e.g. in Figure 2b.\nThe text relies at times heavily on the extended data and figures. Especially the part lines 181 - 216. This makes me wonder\nwhether this manuscript is suited for a short format publication such as the one the authors have chosen to submit to. In this\nsection I am also confused about the CP and EP references. The authors do not separate CP from EP Ninos but rather just\nthe locations over the entire time period. What is the sensitivity of this analysis if it is separated into times of EP and CP\nNinos, and La Nina vs El Nino? Only part of the analysis is split into El Nino and La Nina periods.\nGeneral minor remarks (which are personal suggestions that I believe would increase the readability of the work):\nFor the figure captions please do not say “same as fig x” if the figure is a different one. It is hard to have to go back and forth\nbetween figures to find out what we’re looking at.\nDWC is not a good acronym in my opinion. It would be better to actually name the process since it’s so crucial to this work\nand diabatic cooling / warming has both sides, so the acronym sounds kind of weird to me.\nMinor comments:\nLl 51ff do the authors mean the recharge-discharge mechanism?\nL57 the references are too old, there are more modern studies looking at diabatic processes and how they contribute to El\nNino and La Nina variability\nL 71 the reference to Fig 1a is not very helpful here, the figure just shows the one point where data is available, not the\ndifficulties in sampling.\nL 88 again 5S-5N is mentioned while in the extended data it is stated that the calculations of mixing events are only done\n3S-3N, this needs to be consistent.\nL90 “rapidly diminishes” I think this is an overstatement. It pretty much linearly decreases from the maximum.\nL 93 the reference to barrier layer mixing events either needs to be explained better or not necessary. This speaks to the\noverall a bit “loose” referencing to mixing in this manuscript, this analysis per definition of the method only analysis shear\ndriven mixing, how is this related to the barrier layer mixing events? Related to the “loose” references: How can the shear\ndriven mixing be present at say 4S where there is no strong vertical current shear? The reader needs to be explained what\ntype of mixing we are looking at and why it can arguably be assumed to be active in the region of the study.\nL 95 wind forcing and shear instability might induce different types (at different depths) of mixing. Is deep cycle turbulence\ndistinguished from wind driven turbulence in the upper ocean boundary layer? Are both meant? Is the effect of wind on DCT\nmeant? More precise language would help the reader.\nL 95 - 106 I don’t think I fully follow why this discussion is here. Can you please explain?\nL 111-114 instead of just comparing the strength of the shear Ri should be shown to support this argument. Arguably Ri is\nhard to estimate from mean fields.\nL 115 The explanation of the mixing detection comes too late in the manuscript in my opinion.\nFig 2e CP and EP regions should be indicated\nL187-190 should cite references for the intensification of the EUC during La Nina\nFig 3 shows the occurrence of mixing in 5degree latitude blocks though we have established that the analysis only (at most)\nworks in 3S-3N. What is the sensitivity of the result to reducing the meridional extent? Additionally: why is only the CP\nregion shown?\nFig 4 a I don’t know if I find this figure outstanding, I find the interannual variability quite hard to see. It might help to add an\nindicator of SST anomaly in Nino3.4 or similar to the figure. Over which area are the heat fluxes calculated? Fig 4 now\nanalysis El Nino and La Nina separately, which is nice. I still think that an addition of occurring EP vs CP Ninos is\nnecessary. What if the cooling effect in the east actually represses the before dominant EP Ninos?\nData and Method section\nWhy are GHRSST data used for SST and not ERA5 to be consistent?\nL 470 what is the sorted density?\nJohnson, Gregory C., et al. \"Direct measurements of upper ocean currents and water properties across the tropical Pacific\nduring the 1990s.\" Progress in Oceanography 52.1 (2002): 31-61.\nKarnauskas, Kristopher B., et al. \"The Pacific equatorial undercurrent in three generations of global climate models and\nglider observations.\" Journal of Geophysical Research: Oceans 125.11 (2020): e2020JC016609.", "reviewer_2": "Review of manuscript NCOMMS-24-51953-T\n“Subsurface Ocean Turbulent Mixing Enhances Central Pacific ENSO as Positive Feedback”\nby Liu et al.\nAuthors suggested that importance of turbulent mixing in the equatorial Pacific for CP ENSO, and positive feedback\nassociated with this mixing. I think that authors’ idea is interesting. However, I would like to suggest major revision because\nof following major comments.\nMajor comments:\n1. Authors discuss mixing using Argo float data. Although the number of observations by Argo floats exceeding 44,000 in the\ntropical Pacific after 2000, I do not think that resolution in time and space is enough to discuss turbulent mixing. For\nexample, according to Fig. 1b, the number of observations in the grid of 5 degrees x 3 degrees is around 1000 from 2000 to\n2022. This number means that there are only four profiles per month in the grid of 5 degrees x 3 degrees. Spatial resolution\n(of 3 degrees in meridional) is also too rough to resolve equatorial current system such as Equatorial Under Current and\nTsuchiya jet. I think that these current system, which authors missed, influence mixing around the EL. On the contrary,\nspecial and time scale of vertical mixing is very small such as braking of internal waves, which cannot be resolved by Argo\nfloat data.\n2. Authors should estimate sampling error of Jq, etc, derived from use of Argo data.\nOther comment.\n1. Line 86-87, and Fig.2a. Why occurrence of EL mixing south of the equator is large than that of north of the equator? Is this\nrelated with the sampling number of Argo float profiles?\n2. Caption of Fig.1. Symbols, (cid:1)(cid:2)(cid:3) and (cid:1)(cid:2)(cid:4)(cid:5) , suddenly appeared without explanation."}, "author_response": "Quick response to Reviewer:\nDear Reviewer, thank you very much for the positive evaluation of our manuscript, and\nfor the careful thoughts and comments. The comments and suggestions are helpful for\nus to clarify our results and conclusions. We have carefully studied your comments and\nsuggestions, and revised the manuscript accordingly. We hope you find this version\nmuch clearer than the original one.\nThe responses are written in blue, and the corresponding revisions are made in the\nmanuscript, which are indicated by red line numbers of the new version within the\nresponses.\nMajor comments\nQ1: The authors rely on the state of marginal stability as given by Richardson number\naround 0.25 to substitute stratification for shear in their relationships they employ to\nestimate Jq. The authors rightly state in the extended section that the assumption\nof marginal instability only holds in the presence of strong background currents inducing\nvertical shear. In this manuscript this is assumed to be true within 3°S-3°N around the\nequator. However, the observed undercurrent which gives rise to the strong current\nshear is shown in for example Johnson et al (2002) or Roger (Karnasukas et al 2020) to\nspan at most 2S-2N. I) A rigorous analysis of which profiles, at which location are\nanalyzed needs to be presented, and II) reasonable arguments need to be made why the\nmarginal stability assumption holds, even away from the strongest current system. III)\nWhitt et al (as referenced in the manuscript) examine 3°N partially in relation to strong\nmixing induced by tropical instability waves, no mention of TIW kinetic energy is made\nhere.\nResponse: Thank you for the important questions and suggestions. We first reply to\nyour Point II. We realized that we made some inappropriate statements about the MI\nin the original manuscript.\nFirstly, the assumption of MI (Ri~=0.25) indeed is required, but “strong” vertical\nshear is not necessarily needed; what it really depends on is the relative intensity of the\nshear with respect to the stratification. In other words, if the shear is moderate and the\nstratification is also moderate, MI condition will be fulfilled, too. Nevertheless, in any\ncase, background shear is needed and plays as a fundamental a precondition for shear\ninstability. This should be the correct statement for MI. (We mentioned “In a stratified\nstate, a background vertical velocity shear plays as a necessary condition for MI” in the\noriginal version and keeps it in Line 545 of the new version.)\nSecondly, the MI condition needs not be hold all the time - when this condition is\nfulfilled (maybe occasionally and in a short period), it provides likelihood for\nturbulence to occur, in contrast, when it is not fulfilled, turbulence would less likely to\noccur. The turbulence occurrence would be reflected by the Argo profiles with density\noverturns. Once turbulence occurs, it is assumed the MI condition is met. This could be\nthe cases in the off-equatorial region, where the EUC-induced shear enhancement\ndisappears.\nRegarding the MI in the off-equatorial region, we know that the westward flowing\nSouth Equatorial Current (SEC), which covers a meridional range of about 5°S to 5°N,\nprovides the background shear in the upper ~50 m of the equatorial Pacific region\n(clearly seen in Johnson’ 2002 paper). The SEC is the primary current that provides the\nshear at the MLB and EL depths. Though the shear in 2°S-2°N region is strengthened\nby existence of EUC, the shear in 2°N-3°N and 3°S-2°S are still significant. In the\nmeanwhile, the other factor that needs to be considered is the strength of stratification\nnear the MLB. Also from Johnson’s observations, we see that the stratification is weak\naround the MLB, providing the possibility to form MI condition. The SEC shear could\nbe strengthened or weakened by the zonal wind, which results in higher or lower\nlikelihood of MI and turbulence in the ELs. We show the Ri along 3 meridional sections\nbased on Johnson 2002’s data for your reference:\nIn your Point III you mentioned TIWs, thank you for this suggestion. TIWs are\nenergetic processes in the central to eastern equatorial Pacific. According to our own\nstudies, the Yanai wave mode cover 4°S-4°N associated with strong zonal velocity\ncomponent, and the Rossby wave mode cover 5°S-8°N, associated with strong zonal\nvelocity component particularly at 0°-5°N (Liu et al., 2019) – the zonal velocities\nprovide a significant contribution to the zonal velocity shear in the upper layers (Liu et\nal., 2016); in addition, the Yanai wave mode can also provide weaker stratification in\nparticular phases, providing likelihood of MI from the stratification perspective (Fang\net al., 2023).\nWe added the above arguments to the Data and Method section in the new version\n(Line 638-653).\nThe mentioned references:\nLiu, C., Köhl, A., Liu, Z., Wang, F. & Stammer, D. Deep-reaching thermocline\nmixing in the equatorial pacific cold tongue. Nat Commun 7, (2016).\nLiu, C., Wang, X., Köhl, A., Wang, F. & Liu, Z. The Northeast-Southwest Oscillating\nEquatorial Mode of the Tropical Instability Wave and Its Impact on Equatorial\nMixing. Geophys Res Lett 46, 218–225 (2019).\nFang, L., Ma, K., Liu, C., Wang, X. & Wang, F. Ocean mixing induced by three-\ndimensional structure of equatorial mode tropical instability waves in the Pacific\nOcean. Front. Mar. Sci. 1–14 (2023) doi:10.3389/fmars.2023.1228897.\nNow we come back to your Point I. According to your suggestion and the above\nargument, we checked how many Argo profiles and how frequent turbulence are found\nacross the equator in a greater detail. We first added a panel to Fig. 1b to show the Argo\nprofile numbers in zonal mean at 1° grid in meridional direction. It shows that the Argo\nprofiles are concentrated to the equator; particularly, from the 3°S-2°S band to 2°N-\n3°N band between 170°E-110°W, the accounts of Argo profiles are 3521, 5898, 8259,\n10438, 7116, and 5309, respectively. We conducted the mixing occurrence analysis on\nall the profiles. We then added another panel to Fig. 2 (Fig. 2c) to show the meridional\ndistribution (at 1° grid) of turbulence occurrence rate for the western, central-to-east,\nand far-east equatorial Pacific regions, which shows a meridionally coherent\noccurrence rate in the central-to-eastern Pacific. In our analysis, we didn’t distinguish\nthe meridional variations, instead, we take the 3°S-3°N band as a whole. Whereas, we\nconducted the analysis of mixing coefficients (𝜀 , 𝑘 , 𝐽𝐸𝐿) on all profiles only in\n𝐸𝐿𝑀 𝐸𝐿𝑀 𝑞\nthe 3°S-3°N band where the EL mixing method is valid.\nWe clarified this point in the main text of the new version (Line 87-94 and Line\n164-165).\nQ2: In the introduction of the text the authors reference the number of Argo profiles\n5°S-5°N. These discrepancies need to be reconciled. Why are we shown the number of\nArgo profiles in the 3-5°N/S bands, if they are (rightly) not analyzed?\nResponse:\nDear Reviewer, thank you for this question. We revised the manuscript thoroughly\nbase on this question. Firstly, we keep to show the Argo profiles or the turbulence\noccurrence rate in a wider band in Fig. 1b, Fig. 2a, and Fig. 3a in the new version, to\nindicate the contrast of the turbulence occurrence in the equatorial band to the off-\nequatorial band. When we do this, we intend to emphasize the properties in 3°S-3°N.\nSecondly, in all other places, we focus our analyzation on 3°S-3°N. Some statistics\nthus are reconducted and some panels replotted in the new version, including Fig. 2b,\nc, d and Fig. 3c and f. We note that, the conclusion of the present study is not changed.\nQ3: It is also not very clear in the body of the text how the density overturns relate to\nthe shear driven mixing. This is explained in the extended data, but should be mentioned\nconcisely in the main text as well.\nResponse: Thanks for the suggestion, we added a simple explanation to this point to the\nmain text. Line 87-89.\n“The mixing events are identified as density overturns21–23 , because the latter\nis one of the key features during shear instability-associated billows and turbulences in\nthe shear flow21.”\nQ4: Related: In Figure 3 we are looking at the response of mixing indicators to wind\nstress, which is again not entirely clear to me why. It needs to be better explained in the\ntext how this relates to shear driven turbulence.\nResponse: We also give a more detailed explanation for the wind’s driving effect on the\noccurrence of EL mixing. Line 116-118; Line 230-243.\nQ5: The discussion about the Central Pacific (CP) and Eastern Pacific (EP) is at times\nnot clear. In the abstract it is said that the CP exhibits a stronger signal, what is meant?\nHistorically EP ENSOs are stronger than CP ENSOs. Do the authors mean that CP is\nmore frequent in the past two decades?\nResponse: We are sorry for the confusion. In the revised version, we have made the\ndescription clear. For this specific question, we did mean the frequency of CP ENSO\nin the past two decades. In the new version, we changed it from “while exhibiting more\npronounced signals in the CP than the EP, … in this century” to “while exhibiting\npronounced signals more frequently in the central Pacific than the eastern Pacific in this\ncentury, …”. (Line 26-28).\nQ6: Further, the authors speak about CP and EP but do not show coincidence of mixing\nduring EP or CP events, just the regions.\nResponse: Dear Reviewer, this is a very good question. We agree with you that it is\nbetter to separate EP and CP ENSOs and show whether the mixing’s effect is the same\nor different during two types of ENSO. Unfortunately, we are unable to do so currently\nbecause we don’t have enough pure EP ENSO samples during the Argo era. We\nencountered 7 El Ninos since this century (2002/03, 2006/07, 2009/10, 2014/15,\n2015/16, 2018/19, 2023/24), but only the later 5 El Ninos occurred during the time of\nenough Argo profiles. Moreover, among the 5 El Ninos, none of them belongs to a pure\nEP type, instead, the 2009/10, 2014/15, and 2018/19 El Ninos belong to CP type, while\nthe 2015/16 and 2023/24 El Ninos belong to mixed-type. Therefore, we can’t\ndistinguish CP ENSO from pure EP ENSOs with our data.\nHowever, we noticed that, beside the 3 CP El Ninos, all the other 2 El Ninos\n(including the CP El Ninos）show a CP component, which manifests as cores of positive\nSST anomaly and eastward zonal wind anomalies at the CP region; in the meanwhile,\nthey also show a weak (somehow insignificant) westward zonal wind anomaly at the\nEP region, resembling the 3 CP El Ninos. These features indicate that the pattern of the\ncontrolling factor for the EL mixing, i.e., the wind, of the El Ninos are in general similar\nto a typical CP El Nino. The same happens to the La Ninas: during other pure EP type,\nmixed-type or CP type La Ninas, the wind always shows a westward zonal anomaly in\nthe CP region and weak (insignificant) eastward zonal anomaly in the EP region,\nresembling that of a typical CP type La Nina. Therefore, the present analysis on “CP\nregion” can be regarded as that on a “CP type”.\nFinally, we conclude that, the derived diabatic positive feedback mechanism is\napplicable to CP ENSOs, that’s the reason why we would call it a CP ENSO diabatic\npositive feedback mechanism. In contrary, as discussed above, the diabatic positive\nfeedback should also applicable to the western part of the EP ENSO, yet not to its\neastern part, so that we can’t call the diabatic effect as a positive feedback mechanism\nto the EP ENSO – it becomes only a regional effect to the EP ENSO.\nIn both the text and Figures of the revised version, we pay special attention to the\ndifference between “CP ENSO” and “CP region” to make the description accurate. We\nalso add the above discussion to the Discussion section to justify the proposed idea Line\n354-382.\nWe would like to thank you again for helping clarify our statement.\nQ7: Side note / related minor comment: in some figures it says that CP and EP regions\nare marked, but I cannot see the marking, e.g. in Figure 2b.\nResponse: We added the region indicators such as Nino3, Nino4, and Nino3.4 to this\npanel. After rearrangement, it is Fig. 2f in the new version. In this Section, we would\nlike to only indicate the conventional Nino index regions: Nino3, Nino4, and Nino3.4.\nIn the following more relevant section, 2e specify the CP and EP regions in Fig. 3c and\nf.\nQ8: The text relies at times heavily on the extended data and figures. Especially the part\nlines 181 - 216. This makes me wonder whether this manuscript is suited for a short\nformat publication such as the one the authors have chosen to submit to.\nResponse: Thank you for the comment. Based on both of yours and the Editor’s\nsuggestions, we moved necessary contents and figures particularly for the original\nL181-216 part (Fig. 4e, f, g and h; Fig. 5b and c) from the Extended Data to the main\ntext, so that the text reads smoother now. We also regrouped and rearranged the figure\npanels, which sound more logical and in line with the text. Now the figure numbers in\nthe main text increase from 5 to 7.\nQ9: In this section I am also confused about the CP and EP references. The authors do\nnot separate CP from EP Ninos but rather just the locations over the entire time period.\nWhat is the sensitivity of this analysis if it is separated into times of EP and CP\nNinos, and La Nina vs El Nino? Only part of the analysis is split into El Nino and La Nina\nperiods.\nResponse: It is becoming increasingly clear that this is a good question. We have\nresponded in the response to your Q6.\nQ10: General minor remarks (which are personal suggestions that I believe would\nincrease the readability of the work):\nFor the figure captions please do not say “same as fig x” if the figure is a different one.\nIt is hard to have to go back and forth between figures to find out what we’re looking\nat.\nResponse: This is a good suggestion; we accepted your suggestion and modified the\ncaptions. The revision appears only in the caption of Fig. 3d.\nQ11: DWC is not a good acronym in my opinion. It would be better to actually name\nthe process since it’s so crucial to this work and diabatic cooling / warming has both\nsides, so the acronym sounds kind of weird to me.\nResponse: Thank you for your suggestion, we agree with you and revised accordingly.\nIn the new version, we never use the acronym anymore. Once we need, we use the full\nspell of them, like “diabatic cooling”, “diabatic warming” and “diabatic arming/cooling\nmechanism” directly.\nMinor comments:\nLl 51ff do the authors mean the recharge-discharge mechanism?\nResponse: Yes.\nL57 the references are too old, there are more modern studies looking at diabatic\nprocesses and how they contribute to El Nino and La Nina variability\nResponse: Yes, you are right. In the original version, we only cited studies from the\nhydrologic data analysis. Now in the new version, we added some model-based work\nin recent years. To our best knowledge, our work is the first data-based estimate of the\ndiabatic effect on the basin scale. Our work on the one hand reveals the positive effect\nof turbulent mixing to CP ENSO, on the other hand plays as the gauge for evaluating\nmodel skills. We will be glad to accept your suggested references that are not noticed\nby us in this line. The added references include:\nZhang, R. H., Kleeman, R., Zebiak, S. E., Keenlyside, N., & Raynaud, S. (2005). An\nempirical parameterization of subsurface entrainment temperature for improved\nSST anomaly simulations in an intermediate ocean model. Journal of Climate,\n18(2), 350–371. https://doi.org/10.1175/JCLI-3271.1\nZhang, R. H., Busalacchi, A. J., & DeWitt, D. G. (2008). The roles of atmospheric\nstochastic forcing (SF) and oceanic entrainment temperature (Te) in decadal\nmodulation of ENSO. Journal of Climate, 21(4), 674–704.\nhttps://doi.org/10.1175/2007JCLI1665.1\nZhang, R. H., Zheng, F., Zhu, J., & Wang, Z. (2013). A successful real-time forecast\nof the 2010-11 la Niña event. Scientific Reports, 3(January).\nhttps://doi.org/10.1038/srep01108\nGuan, C., Tian, F., McPhaden, M. J., Hu, S. & Wang, F. Zonal Structure of Tropical\nPacific Surface Salinity Anomalies Affects the Eastern and Central Pacific El\nNiños Differently. Geophys Res Lett 50, (2023).\nL 71 the reference to Fig 1a is not very helpful here, the figure just shows the one point\nwhere data is available, not the difficulties in sampling.\nResponse: Yes, thank you, we deleted the reference here.\nL 88 again 5°S-5°N is mentioned while in the extended data it is stated that the\ncalculations of mixing events are only done 3°S-3°N, this needs to be consistent.\nResponse: Thanks for the suggestion, we agree with you that the mentioned quantity\nindexes should be consistent regarding the region. Accordingly, we made revisions in\nthe new version, please see our response to your Q2.\nL90 “rapidly diminishes” I think this is an overstatement. It pretty much linearly\ndecreases from the maximum.\nResponse: We delete “rapidly” from the text.\nL 93 the reference to barrier layer mixing events either needs to be explained better or\nnot necessary. This speaks to the overall a bit “loose” referencing to mixing in this\nmanuscript, this analysis per definition of the method only analysis shear driven mixing,\nhow is this related to the barrier layer mixing events? Related to the “loose” references:\nHow can the shear driven mixing be present at say 4S where there is no strong vertical\ncurrent shear? The reader needs to be explained what type of mixing we are looking at\nand why it can arguably be assumed to be active in the region of the study.\nResponse: Sorry for this loose reference. We have deleted this reference and this\nsentence. [We will analyze the similarity and difference between the two patterns in the\nfuture.]\nL 95 wind forcing and shear instability might induce different types (at different depths)\nof mixing. Is deep cycle turbulence distinguished from wind driven turbulence in the\nupper ocean boundary layer? Are both meant? Is the effect of wind on DCT meant?\nMore precise language would help the reader.\nResponse: Thank you for the question and suggestion. The description of the two\nmechanisms for the generation of the EL mixing is greatly expanded. Line 116-148.\nL 95 - 106 I don’t think I fully follow why this discussion is here. Can you please\nexplain?\nResponse: Sorry for the confusion. This paragraph, by showing the correlations\nbetween the two variables, meant to explain the wind-driven mechanism for the\ngeneration of the EL mixing. In the meanwhile, the following paragraph explains the\nshear instability mechanism for the generation of the EL mixing. In the new version,\nwe rearranged the paragraphs (see the above Response), so that this discussion becomes\nclearer than the original version. Line 116-148.\nL 111-114 instead of just comparing the strength of the shear Ri should be shown to\nsupport this argument. Arguably Ri is hard to estimate from mean fields.\nResponse: Thank you for this suggestion. We followed your suggestion and plotted the\nRi that is calculated from the mean S2 and mean N2, which confirms our argument. We\nadded the mean S2, N2 and Ri panels to Extended Data Fig 1, and revised the main text.\nLine 132-148.\nL 115 The explanation of the mixing detection comes too late in the manuscript in my\nopinion.\nResponse: We understand your feeling. Accordingly, in the new version, we move 2\nsentences of the mixing detection method to the main text, but would keep most content\nof the mixing detection method in the Data and Method section to keep the main text\nconcise. Line 87-97.\nFig 2e CP and EP regions should be indicated.\nResponse: In this Section, we only indicate the conventional Nino index regions: Nino3,\nNino4, and Nino3.4. We specify the diabatic effect-determined CP and EP regions in\nthe late figures. In the revision, we clarified that the specific CP (170°E-140°W) and\nEP (130°W-110°W) regions as defined in our study is based on the diabatic effect, we\nalso compared them with the conventional Nino4 and Nino3 regions in terms of both\nthe area and the meaning behind them (Line 218-229). Thank you for your comments\nthat lead to the clarification of the definitions of the CP and EP regions.\nL187-190 should cite references for the intensification of the EUC during La Nina.\nResponse: Thank you for the suggestion. The argument was based on results of our own\ncalculation with TAO data at 170°W, 0°N; the related figures are displayed as Extended\nFig. 5c-f. In the original version, we cite the corresponding panels too late. In the\nrevised version, we cite the panels timely to avoid confusion. We also cite a reference\nwith the same conclusion based on 140°W, 0°N (Warner and Moum, 2019).\nFig 3 shows the occurrence of mixing in 5° latitude blocks though we have established\nthat the analysis only (at most) works in 3°S-3°N. What is the sensitivity of the result\nto reducing the meridional extent? Additionally: why is only the CP region shown?\nResponse:\n1) In the revised version, according to your questions here and in the Major Comments,\nwe use only quantities in the 3°S-3°N region to do the analysis, and find that the\nconclusion is not altered. Please see our response to you Q2.\n2) That of the EP region was also plotted but was put in the Extended Data Figure 6 in\nthe original version. In the revised version, we put the analysis of both CP and EP\nregions on a same figure, the new Figure 4.\nFig 4a I don’t know if I find this figure outstanding, I find the interannual variability\nquite hard to see. It might help to add an indicator of SST anomaly in Nino3.4 or similar\nto the figure. Over which area are the heat fluxes calculated?\nResponse: We replotted this panel and rearrange the figures. In the new version, the Fig.\n5 and Extended Data Fig. 3 show the interannual variations. In the meanwhile, we show\nthe ONI index aside both figures. This makes the depiction and description clearer than\nthe original version. Thank you for the comment.\nFig 4 now analysis El Nino and La Nina separately, which is nice. I still think that an\naddition of occurring EP vs CP Ninos is necessary.\nResponse: Please see our responses to your Q6.\nWhat if the cooling effect in the east actually represses the before dominant EP Ninos?\nResponse: Thank you for this question. In my opinion, the repression (mostly from the\nsurface heat flux’s damping effect) would not prevent the growth of an EP El Nino,\nbecause the dominant and active mechanisms for EP El Nino’s growth is the wind-\ndriven thermocline mechanism and zonal advection mechanism in the EP region, while\nthe surface heat flux’s damping effect is a passive mechanism which reflects the SST\nchange. In the meanwhile, the positive diabatic CP ENSO dynamics would help with\nthe growth of the western part of the EP EL Nino.\nData and Method section\nWhy are GHRSST data used for SST and not ERA5 to be consistent?\nResponse: We used GHRSST to represent observation. In fact, they don’t have obvious\ndifferences from ERA5 at the interannual time scale.\nL 470 what is the sorted density?\nResponse: Due to EL mixing, the density manifests as density overturns, which means\nthe density is not monotonically increasing with depth. In this case, firstly the density\nis sorted, and then N2 is calculated based on the sorted density.\nThank you again for your careful and instructive comments, which help improve\nour manuscript.\nJohnson, Gregory C., et al. \"Direct measurements of upper ocean currents and water\nproperties across the tropical Pacific during the 1990s.\" Progress in Oceanography 52.1\n(2002): 31-61.\nKarnauskas, Kristopher B., et al. \"The Pacific equatorial undercurrent in three\ngenerations of global climate models and glider observations.\" Journal of Geophysical\nResearch: Oceans 125.11 (2020): e2020JC016609."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I have reviewed the changes by the authors and am almost satisfied. I commend the authors for a truly interesting and\ninspiring study. One can always argue that there is too little data to draw conclusions, but I believe this examination of the\navailable data is very noteworthy and deserves to be published in Nature Communication. I believe the manuscript has\nimproved and is almost ready for publication. The one important point I am not satisfied with is the response to the distinction\nbetween wind and \"shear\" /MI (Ri, the authors point out that shear is only one component to get to Ri=0.25 which is correct,\nbut in a stratified upper ocean it is the main one that changes and therefore moves Ri from above to below 0.25) driven\nmixing. I quote the author's rebuttal:\nQ4: Related: In Figure 3 we are looking at the response of mixing indicators to wind stress, which is again not entirely clear to\nme why. It needs to be better explained in the text how this relates to shear driven turbulence.", "reviewer_2": "Because authors well responded to my comments, I suggest that this paper is accepted in Nature Communications after revision for following tiny comments.\n1. Line 118. What is “MIs”?\n2. Line 918. “feedback oo ENSO” -> “feedback to ENSO”?"}, "author_response": "Dear Reviewer:\nThank you very much for further thinking this issue. It is a good suggestion to distinct the role\nof wind stress from the role of strong shear/MI in forming the turbulent mixing around the mixed\nlayer base. We here try our best to answer your question in the following 3 subsections, which are\nfollowed by the revision to the text according to your suggestion.\n1, Two mechanisms for wind stress in promoting the generation of turbulent mixing around the\nmixed layer base:\n1) Mechanism I: providing a background shear instability condition. The westward zonal wind\nstress (of the trade wind) maintains the westward south equatorial current (SEC) in the near-\nsurface layer via momentum input and the eastward equatorial undercurrent (EUC) in the\nsubsurface layer via west-east pressure (Yu and Mcphaden, 1999), which jointly form a\nbackground vertical shear above the EUC core. The shear of (primarily) SEC (and possibly also\nEUC) beneath the mixed layer base provides the basic precondition for shear instability (random\nvorticity and/or internal wave perturbations then may initiate the turbulence). In general, the\nstronger the trade wind stress, the stronger equatorial currents and the resulted vertical shear,\nand the higher the possibility of turbulence generation. Therefore, it is obviously seen that the\nvariations of zonal wind stress at longer (>monthly) time scales lead to the variations of both\nthe shear (Supplementary Fig. 5abcd) and the mixing indicator (Fig. 2d, 3cf, 4bf) at the same\ntime scales. This mechanism was explained in the Line 230-243 in the last version.\n2) Mechanism II: wind-driven entrainment. The wind stress provides momentum input to the\nmixed layer, breaks the equilibrium of the ML and leads to entrainment at the mixed layer\nbottom, as well shown in both laboratory experiments (e.g., Strang and Fernando 2001) and\ndirect numerical simulation (DNS) studies (e.g., Pham and Sharkar, 2017). Entrainment is\nfeatured as a process of wind-driven turbulence originated at the mixed layer bottom deepening\nthe mixed layer into the depths of the original stratified shear flow.\n2, Two mechanisms for strong shear/MI in the generation of turbulent mixing beneath the mixed\nlayer base.\n1) Strong shear or MI or low Ri below the mixed layer base favors the occurrence of shear\ninstability and turbulence, which may be initiated by suitable perturbations from any\nmechanisms.\n2) This condition also favors the entrainment to develop and to penetrate to a deeper layer, forming\nthe observed DCL in the central-to-eastern equatorial Pacific.\n3, Interwind roles of them under MI conditions, and separated roles in specific conditions\n1) It should be clarified that the roles of wind stress and shear are not independent, particularly\nunder the MI condition which is frequently fulfilled in the central equatorial Pacific (Moum et\nal., 2023). This is what our study involves, it is also what many previous works was focusing\non. On the one hand, wind stress controls the strength of shear, which may easily move a\nmarginal instability condition to a stability or instability condition; on the other hand, strong\nshear favors the wind-induced entrainment to develop into DC turbulence, and favors the\nentrainment-associated internal wave into shear instability to occur as well, and vice versa. As\nsuch, the MI provides the precondition for the occurrence of turbulence, while the wind stress\ninitiates the turbulence there– via providing additional shear, its induced internal wave and/or\nentrainment at the mixed layer bottom. Moum et al. (2023) observed that stronger wind stress\nwill lead to larger and quicker enhancement of mixing beneath the mixed layer base. This is the\nreason why the proposed parameterization for mixing coefficients constitutes both wind stress\nand shear.\n2) In contrast, imagine two scenarios where the wind stress and the MI condition are separated\nfrom each other. One is a strong thermocline with weak shear or large Ri just beneath the mixed\nlayer base, where additional wind stress may lead to entrainment that will be confined in a very\nthin layer and can’t support a DC turbulence. If this is the case in the study area, the entrainment-\nassociated turbulence patch might be very small and can’t be detected by Argo profile data. The\nother is a low Ri/MI condition far away from the surface layer, like the observed upper core\nlayer (UCL, Moum et al., 2009; Inoue et al., 2012), where the shear instability and turbulence\nis found by ourselves generated by local dynamics and has nothing to do with wind stress (Liu\net al., 2020). If this is the case in the study area, it is excluded because we focus only on mixing\naround the mixed layer base.\n4, For the revision of the text:\nWe revise this section by providing the individual mechanisms of them to distinct their roles in the\ngeneration of EL mixing. We also discuss the complicated interactions between wind- and shear-\ninduced mixing by providing an additional paragraph.\nWe modify the wind stress’s role with a long sentence:\n“Firstly, wind stress forcing offers one of the basic interpretations for occurrence\nof EL mixing. Stronger zonal wind stress (stronger trade wind) leads to higher\noccurrence of EL mixing via two mechanisms. The first is it provides stronger vertical\nshear and hence more shear instability condition in subsurface layer, which will be\ndemonstrated in the following subsection, and the second is it leads to stronger and\nquicker enhancement of entrainment turbulence, which are demonstrated by laboratory\nexperiment (Strang and Fernando 2001), direct numerical simulations (Pham and Sharkar,\n2017), and observed by several years of turbulence measurements (Moum et al., 2023).”\nSimilarly, we modify the shear’s role with a longer sentence:\n“Secondly, local background shear stability condition20 offers another basic\ninterpretation for occurrence of EL mixing. The occurrence of EL mixing reflects\nfulfillment of local shear instability therein.”\nMoreover, we also added a paragraph below those two paragraphs to clarify their combined roles:\n“It is noted that the roles of wind stress and shear in promoting the generation of\nEL mixing are not independent, particularly under the MI condition which is frequently\nfulfilled in the central-to-eastern equatorial Pacific region (170°E–110°W). On the one\nhand, wind stress controls the strength of subsurface background shear and the\ninitiation of entrainment at the MLB, on the other hand, strong background shear or\nMI favors shear instability and the wind-driven entrainment to develop into deep cycle\nturbulence. Overall, the MI in the subsurface layer provides the precondition for the\noccurrence of turbulence, while the wind stress initiates the turbulence there.”\nWe hope we have addressed your concern and made the manuscript clearer. Thank you again.\nReferences in the response:\nLiu, C., X. Wang, Z.Liu,A.K¨ohl, W. D. Smyth, and F. Wang, 2020 (this is the leading author’s\nwork): On the formation of a subsurface weakly sheared laminar layer and an upper thermocline\nstrongly sheared turbulent layer in the eastern equatorial Pacific: Interplays of multiple-time-scale\nequatorial waves. J. Phys. Oceanogr., 50, 2907–2930.\nMoum, J. N., R. C. Lien, A. Perlin, J. D. Nash, M. C. Gregg, and P. J. Wiles (2009), Sea surface\ncooling at the equator by subsurface mixing in tropical instability waves, Nat. Geosci., 2, 761–765.\nInoue, R., R.-C. Lien, and J. N. Moum (2012), Modulation of equatorial turbulence by a tropical\ninstability wave, J. Geophys. Res., 117.\nPham, H. T. and S. Sarkar (2017): Turbulent entrainment in a strongly stratified barrier layer, J.\nGeophys. Res. Oceans, 122, 5075–5087.\nStrang, E. J., and H. J. S. Fernando, 2001: Entrainment and mixing in stratified shear flows. J. Fluid\nMech., 428, 349–386.\nYu, X., and McPhaden, M. J. (1999). Dynamical Analysis of Seasonal and Interannual Variability\nin the Equatorial Pacific. J. Phys. Oceanogr. 29 (9), 2350–2369.\nReviewer #2:\nBecause authors well responded to my comments, I suggest that this paper is accepted in Nature\nCommunications after revision for following tiny comments.\n1. Line 118. What is “MIs”?\n2. Line 918. “feedback oo ENSO” -> “feedback to ENSO”?\nResponse:\nAll are revise, thanks."}]

Decision letter