The role of upper-ocean variations of the Kuroshio-Oyashio Extension in seasonal-to-decadal air-sea heat ﬂ ux variability

The Kuroshio-Oyashio Extension (KOE) is the North Paci ﬁ c oceanic frontal zone where air-sea heat and moisture exchanges allow strong communication between the ocean and atmosphere. Using satellite observations and reanalysis datasets, we show that the KOE surface heat ﬂ ux variations are very closely linked to Kuroshio Extension (KE) sea surface height (SSH) variability on both seasonal and decadal time scales. We investigate seasonal oceanic and atmospheric anomalies associated with anomalous KE upper ocean temperature, as re ﬂ ected in SSH anomalies (SSHa). We show that the ocean-induced seasonal changes in air-sea coupled processes, which are accompanied by KE upper-ocean temperature anomalies, lead to signi ﬁ cant ocean-to-atmosphere heat transfer during November-December-January (i.e., NDJ). This anomalous NDJ KOE upward heat transfer has recently grown stronger in the observational record, which also appears to be associated with the enhanced KE decadal variability. Highlighting the role of KOE heat ﬂ uxes as a communicator between the upper-ocean and the overlying atmosphere, our ﬁ ndings suggest that NDJ KOE heat ﬂ ux variations could be a useful North Paci ﬁ c climate indicator.

The KE region (31°−36°N, 140°−165°E) is especially characterized by interannual and decadal timescale and large-scale air-sea coupled processes in the midlatitude North Pacific, such as windinduced oceanic Rossby wave propagations 5,12,19,20 . The recent KE modulations since 1980 between the positive and negative phase of KE SSH anomalies have preferred decadal time scales 8,15,[21][22][23] . To be specific, a positive KE phase (e.g., positive KE SSH anomalies, Fig. 1a) is accompanied by an enhanced southern recirculation gyre, stronger KE eastward transport, and a northward-shifted oceanic jet 9,15 . The detailed description and synthesized schematic can be found in Fig. 1 of Qiu (2019) 24 . During the positive KE phase, the intensified KE current and eastward jet enhance tropical-origin warm water transport into the KE region and beyond into the midlatitude North Pacific, thereby significantly impacting ocean thermodynamic properties 5,9 . The reverse holds when the KE transitions to the negative KE SSH phase.
The decadal KE SSH variability is known to be associated with large-scale Pacific climate modes -constituting the Pacific Decadal Oscillation (PDO) and North Pacific Gyre Oscillation 6,9,14,15,22,25,26 . In addition, observational and modeling studies have suggested that midlatitude KE SSH fluctuations may closely interact with the subtropical and tropical variability like the Pacific Meridional Modes (PMM) and El Niño-Southern Oscillation (ENSO) 21,27 . Those studies highlight the KE downstream response as a slow nudging of atmospheric circulation anomalies that can reduce/enhance the trade winds, which are key to driving changes in the subtropical and tropical Pacific variability. Also, other studies have described the KE low-frequency expression as a new mode of quasi-decadal Pacific climate variability named ocean-atmosphere coupled Pacific Decadal Procession 28,29 . Moreover, a recent study suggests that the North Pacific SSH anomalies, especially along the KE region, may contribute to initiating ocean temperature extremes like Northeast Pacific Marine Heatwaves 30 . Growing observational evidence suggests that the KE SSH variability could be a crucial climate indicator of extratropical-tropical teleconnections, air-sea coupling, and extreme events, potentially leading to Pacific climate predictability on decadal time scales. However, a fundamental question remains to be clarifiedwhat are the mechanisms and characteristics by which interannual to decadal scale variations of SSH or accompanied anomalies physically impact air-sea fluxes, and how have the seasonal impacts of those strong decadal fluctuations affected the Pacific climate?
In this study, we use KE SSH variability as an indicator of changes in upper-ocean heat content over the Northwestern Pacific 6,8,18,20,[31][32][33][34][35] ; many previous studies have shown that on seasonal-to-decadal scales, changes in upper-ocean heat content account for most of the large-scale changes in SSH 32,33 . As ocean heat content variability is largely governed by temperature over the KE region 6 , we focus on the role of the anomalous upperocean temperatures associated with KE SSH variability in driving regional and large-scale air-sea interactions on seasonal-todecadal time scales (e.g., rather than describing small-scale oceanic processes 36,37 ). We first identify key patterns of surface latent and sensible heat flux anomalies associated with KE SSH anomalies (SSHa). We then investigate seasonal variations in coupled ocean-atmosphere processes, including air-sea heat transfer, due to impacts of upper-ocean temperature anomalies on SST and atmospheric circulation and, finally, show how the earlier winter KOE surface heat flux anomalies are closely linked to KE SSH variability on seasonal and decadal time scales.

Decadal KE upper ocean variations
The cross-correlation between the KE SSH and 0-700 m averaged KE ocean temperature (Fig. 1b) shows that KE SSH anomalies are strongly accompanied by upper-ocean thermal anomalies on decadal time scales 33,38 . These strong decadal variations of subsurface KOE are known to be modulated by basin-scale wind stress curl variability via baroclinic oceanic Rossby wave propagation that induces anomalous thermocline depth and ocean heat content 5,19,31,[39][40][41][42] . The sinusoidal shape of the cross-correlation indicates that those oceanic variations have a decadal oscillatory nature, which corresponds to a period of~9-11 years. Linear regression of zonally-averaged (140°−165°E) ocean temperature anomalies onto the KE SSH index (Fig. 1c) shows the largest temperature anomalies at the KE band (31°−36°N), particularly at 300-600 m depth which is near the mean thermocline. Given the significant contributions of temperature anomalies (e.g., larger contributions than salinity) on the ocean heat content variability over the North Pacific WBC region 6 , we interpret the ocean temperature anomalies associated with KE SSH in this study as the upper-ocean heat content accompanied by wind-induced KE SSH changes. Note that in the rest of the paper, we use the normalized KE SSH index, which is divided by one standard deviation, to show typical magnitudes of KE SSH anomalies associated with upper ocean temperature variations. The KOE region is characterized by a strong SST gradient and isotherms that slopes steeply downward toward the southeast (pink lines in Fig. 1c). Thus, as the latitude of the KOE fluctuates, this region experiences significant seasonal-to-decadal variations in SST, mixed layer depths, and the upper-ocean heat budget 32,43,44 .
We continue to show and compare surface anomaly patterns associated with KE SSH using regression maps of SSH, depthaveraged ocean temperature, SST, and the sum of the latent surface heat flux (LHF) and sensible surface heat flux (SHF) with respect to the KE SSH index ( Fig. 1a and d-f). Those regression maps (with no lags) exhibit that the ocean thermal anomalies (Fig. 1d) show a high correspondence with the SSH anomaly pattern (Fig. 1a) and also with SST ( Fig. 1e) and heat flux (Fig. 1f) anomalies. The same sign of ocean temperature, SST, and heat fluxes (upward positive) anomalies, especially over the eastward KE jet, reveals the export of heat from the ocean to the atmosphere, suggesting that those ocean temperature anomalies did not originate from local atmospheric forcing. We note that while anomalous SSH (Fig. 1a) and ocean temperature (Fig. 1d) are characterized by quasi-stationary meander patterns over the 144°E and 150°E 9 , the SST regression (Fig. 1e) shows a rather complex pattern, including a frontal boundary signature (e.g., North-south dipole anomaly pattern) as well as the local KOE anomalies. We later discuss the role of those large-scale dipole SST anomalies resembling the North-South phase pattern (see Supplementary  Fig. 1) in seasonal KOE evolutions 22 , but here we focus on the local KOE anomalies that represent the anomalous oceanic heat content. The surface heat flux regression (Fig. 1f), with positive anomalies over the upstream KE jet and its downstream, indicates that the upward heat fluxes act to damp the warm SST anomalies.
It is important to note that the most prominent response of KOE LHF and SHF anomalies to KE SSH is found when the KE SSH leads heat flux anomalies by 2-4 months, not at lag 0 months (Fig. 2). Given the delayed response of KOE heat fluxes to KE SSH, we focus on the lagged KOE heat flux anomalies associated with KE SSH and use their spatial structures to identify corresponding temporal evolutions in Fig. 3. In Fig. 3a, the 3-month lagged LHF anomalies regressed on the KE SSH index reveal upward heat flux over the upstream region [140°−153°E], especially the two quasi-stationary meanders, and over the larger downstream region [153°−165°E]. A cross-correlation function between the KE SSH and KOE LHF timeseries (see methods for details) confirms their strong covariability on decadal time scales (Fig. 3b). The oscillatory character and time scale of KOE LHF and KE SSH are nearly identical to those of KOE oceanic variations shown in Fig. 1b (Note that no filtering is applied to the indices). When a 1-year running mean is applied to KOE LHF, the correlation between KE SSH and low-pass filtered KOE LHF reaches 0.7 (Fig. 3c), revealing that near half of the low-frequency KOE LHF variability can be explained by upper-ocean dynamics (e.g., KE SSH or KOE ocean temperature anomalies). Analysis of KOE SHF (Fig. 3d-f) shows similar results of LHF ( Fig. 3a-c), except that the prominent SHF anomalies are found in the relatively northern KOE region poleward of 42 o N (Fig. 3d).

Seasonal KE variations and associated air-sea coupled processes
To investigate seasonal variations of KE upper-ocean variability, we first examine seasonal changes in ocean temperature anomalies associated with KE SSH and their impacts on air-sea coupled  b Cross-correlation function of the monthly KE SSH and KOE LHF indices (no filtering applied). The dashed and solid red lines indicate correlation coefficients statistically significant at the 5% and 10%, respectively, based on a two-sided t-test with the adjusted degrees of freedom. c Monthly time series of the KE SSH index (blue) and 1-year running mean KOE LHF index (red) with a correlation between the two. d-f Same as in (a-c) but for the surface SHF. Note that the SHF amplitude is nearly half of the LHF amplitude (compare a and d). In a and d, the black contours indicate regression coefficients statistically significant at the 10% level, and the black box is the domain used for defining the KOE LHF/SHF indices (see methods for details). processes (Fig. 4). Seasonal regression maps between 3-month running mean ocean temperature and SST anomalies with the 3-month running mean KE SSH index are shown in Fig. 4a, b (e.g., the panel (JFM) indicates the JFM anomalies regressed on the JFM KE SSH index). We note that the upper (0-100 m) ocean temperature anomalies over the KE region become stronger from late spring to summer (AMJ to JJA in Fig. 4a), accompanied by warm SST anomalies over the KOE region (AMJ to JJA in Fig. 4b). The significant variations of ocean anomalies during the warm season can be supported by previous findings that the Kuroshio transport across the south of Japan (e.g., Tokara Strait) and at the upstream region has maximum in spring/summer 24 . As the contribution of advection by the Kuroshio system becomes increasingly important in summer 45 , the heat transport via horizontal advection seems to be a major contributor in affecting upper ocean temperature anomalies associated with the KE SSH. The stronger oceanic connections between the surface and subsurface in spring/summer despite increasing vertical stratification, however, are questionable and will be revisited in future studies. From Fig. 4a, b, we find that the enhanced upper ocean warm anomalies substantially increase the local KOE SST warming from AMJ to JJA, which may affect the overlying atmosphere 46,47 .
We next show the impacts of the increase in summer KOE surface warming on anomalous atmospheric circulations (e.g., geopotential heights) during the KE SSH events (Fig. 4b, c). Seasonal changes in troposphere vertical motions associated with KE SSH variability are shown in Fig. 4c. We note that horizontal and vertical atmospheric anomalies associated with KE SSH show substantial changes from the spring to summer. Specifically, the barotropic vertical structure of height pressure anomalies with a north-south orientation, which was maintained in the late winter and early spring (JFM to MAM in Fig. 4b, c) and accompanied by a basin-scale SST gyre pattern (shadings of JFM to MAM in Fig. 4b), transitions into a baroclinic vertical structure with dipole anomalies tilting poleward after spring (AMJ to SON in Fig. 4c) as the local KOE SST warming becomes prominent (shadings of AMJ to JJA in Fig. 4c). We attribute the changes in geopotential height anomalies to an increase in baroclinic instability arising from the intensified cross-frontal SST gradient due to KOE surface warming, as previous studies have demonstrated [46][47][48][49][50][51][52] . Studies have revealed that nonlinear dynamics, such as transient eddy vorticity fluxes, can drive an atmospheric response to extratropical SST anomalies. In particular, Smirnov et al. (2015) show that the strong upward velocity anomalies of the vertical atmospheric circulation over the warm SST anomalies enhance low-level eddy moisture fluxes and reduce the low-level stability. The increasing lower tropospheric baroclinicity shown in AMJ-SON of Fig. 4c is supported by seasonal changes in the Eady growth rate Fig. 4 Seasonal impacts of upper ocean variability on oceanic and atmospheric anomalies. a Depth-latitude maps of the regression coefficient for the linear regression of zonal-mean (153°-165°E) ocean temperature anomalies (°C) with the normalized KE SSH index, which is divided by one standard deviation, without lags. b Spatial maps of regression coefficients for the linear regression of SST (shading;°C) and 500hPa geopotential height (contours; m) anomalies on the normalized KE SSH index without lags. c Latitude-pressure maps of correlation coefficients between the zonal-mean (140°E-180°E) geopotential height anomalies (m) and the KE SSH index without lags. The black contours indicate regression coefficients statistically significant at the 10% level.
( Supplementary Fig. 2), where the Northwestern Pacific is characterized by strong dipole anomalies from summer to autumn. It is important to note that the increasing low-level baroclinicity appears to lead the subpolar (i.e., northern) lowpressure anomalies to be displaced toward the south of the latitude 38 o N, contributing to a deepening anomalous low in the lower troposphere in early winter (OND and NDJ in Fig. 4b, c). Those two cells-which in our study form a tropospheric circulation with downward (southward) flow over the south of the KOE front and upward (northward) flow over the north of the front (AMJ-SON of Fig. 4c)-are consistent with the atmospheric response to warm KOE SST anomalies reported in Smirnov et al. 2015 (see their Fig. 4c), which has been only detected in a highresolution model. While the detailed development of those largeamplitude low-pressure anomalies during earlier winter (OND and NDJ in Fig. 4b) needs to be further diagnosed, our results show that the strong low-pressure system of OND and NDJ might be closely associated with ocean-driven changes (e.g., KOE upperocean temperature anomalies and related SST warming).
In this study, we suggest that the earlier winter atmospheric anomalies induced by upper-ocean KE variability play a key role in driving significant seasonal KOE air-sea heat exchanges. Consistent with Fig. 4b, Fig. 5a shows noticeable differences in the sea level pressure (SLP) patterns associated with KE SSH between the early (e.g., SON to NDJ) and late (e.g., DJF to MAM) cold seasons. From DJF to MAM, the KE SSH variability is accompanied by a basin-scale dipole SLP pattern, which resembles the North Pacific Oscillation 53 . From SON to NDJ, on the other hand, the northern lobe of low SLP anomalies becomes deeper and shifts towards the south near the KOE downstream (as also shown in Fig. 4b), resembling the Aleutian low 54 . Note that the anomalous upward heat fluxes become prominent during only the early cold season (OND, NDJ, and DJF in Fig. 5a), whereas no significant anomalies are found in any other seasons. We suggest that the 1-month preceding Aleutian-low like anomalies (SON, OND, and NDJ in Fig. 5a), which are associated with increasing baroclinicity due to warming KOE surface in Fig. 4b, c, are important to lead the significant KOE upward surface heat fluxes, via northwesterly winds that bring colder and drier air masses into the Northwestern Pacific.
To measure the ocean's physical contribution to surface heat flux anomalies, we examine the relative contributions of the anomalous wind and sea-air temperature difference (SST-SAT) to heat flux changes (Fig. 5b, c). The anomalous SHF and LHF in Fig. 5 and Supplementary Fig. 3 are estimated as bulk surface fluxes as follows.
In Eq. (1), primes indicate anomalies, ρ is the air density, c p is the specific heat at constant pressure, C s is the sensible heat bulk transfer coefficient, W is the scalar wind speed, 4T is the difference of sea-air potential temperature, L is the latent heat of water vaporization, C L is the latent heat bulk transfer coefficient, and 4Q is the difference of specific humidity. The areal-mean SHF anomalies in the KOE region are nearly identical to the KOE SHF index of this study (R = 0.91 for monthly index), and so we use wind speed and SST-SAT anomalies averaged over the KOE region  Fig. 5b) than between SHF and wind speed (bottom panels in Fig. 5b) reveal that the ocean-atmosphere temperature difference dominantly controls the KOE air-sea heat exchange. The regression map of the NDJ SST (shading in Fig. 5c) and SST-SAT (contour in Fig. 5c) onto the ASO KE SSH index shows that the seaair temperature differences are primarily characterized by the SST anomaly pattern rather than SAT anomalies ( Supplementary Fig. 4), describing the SST as the critical contributor. Our results are consistent with recent findings using a high-resolution ocean general circulation model 45 , where the ocean-atmosphere temperature difference is a crucial predictor for winter variations of the net surface heat fluxes, especially the upward turbulent fluxes. The same analysis but with LHF fluxes (Supplementary Fig. 3) reveals that the SST, by controlling the saturation specific humidity at the ocean surface, also plays a more important role than the wind speed in determining the earlier winter LHF (especially in OND).
Finally, seasonal lead-lag correlations between the KE SSH and KOE heat fluxes (Fig. 6a, b) show that the KE SSH significantly leads the KOE LHF and SHF in satellite records, where the earlier cold season (e.g., OND to DJF) KOE heat fluxes are consistently led by warm season KE SSH. A lead-lag correlation of the sum of KOE SHF and LHF (Fig. 6c) shows their maximum correlations when the summer (especially ASO) KE SSH leads the earlier winter (especially NDJ) KOE heat fluxes. The substantial resemblance between the ASO KE SSH and NDJ KOE heat flux time series (Fig. 6d-f) reveals that the KOE heat flux variations effectively represent the KE SSH variability not only on decadal (Fig. 6) but also on seasonal scales (Fig. 3), explaining a substantial portion (~50%) of its total variance. Specifically, our findings show that NDJ KOE heat fluxes are critically linked to KE SSH fluctuations, which could be a key player linking the KE SSH and other Pacific climate variability that have been observed and discussed in previous studies (e.g., Pacific Decadal Precession 28,29 , ENSO 21,27 , and Northeast Pacific Marine Heatwaves 30 ).

DISCUSSION
Using observational satellite and reanalysis, we have examined the role of upper-ocean variations over the Kuroshio-Oyashio Extension (KOE) region in the Northwest Pacific. While previous literature has extensively focused on the low-frequency component of the KOE oceanic variations (e.g., SSH, ocean heat content, or upper-ocean temperature), the present study centers on seasonal changes in coupled ocean-atmosphere processes of those ocean decadal fluctuations. Consistent with the previous view, our findings support the view that the surface heat fluxes are key indicators of KOE decadal variations 19,20 . We further show that earlier winter (i.e., NDJ) KOE latent (LHF) and sensible (SHF) heat fluxes are the main conveyor of decadal ocean variability to seasonal air-sea heat fluxes.
Our main findings indicate that seasonal ocean heat content anomalies associated with KE SSH variations can be a significant source of earlier winter air-sea heat exchange variations over the KOE. From spring to summer, the KOE upper ocean warm temperature anomalies associated with KE SSH become stronger. By increasing ocean surface warming and the meridional SST gradient over the Northwestern Pacific, the summer ocean temperature anomalies associated with KE SSH appear to enhance the low-level atmospheric baroclinicity, contributing to a significant transition in coupled air-sea structures. Our results, where the North-South dipole SLP pattern evolves into the large-scale monopole atmospheric circulation (e.g., Aleutian-low like pattern) from spring to fall, suggest that the anomalous ocean heat content accumulated during summer might play a role in upward KOE heat transfer during earlier winter. During summer, the relatively weaker wind speeds (i.e., weaker prevailing Northwesterly winds) and the shallower mixed layer depth might allow the subsurface ocean to accumulate contributions of anomalous ocean heat content driven by geostrophic advection (e.g., the arrival of westward propagating Rossby waves). This anomalous ocean heat content accumulated during summer might be especially important for the upward KOE heat flux anomalies during earlier winter (i.e., October-to-December), where the background wind changes are still weaker compared to the late winter (i.e., December-to-February).  anomalies over the KOE 17,55 , our findings indicate that there is a strong consistency between seasonal (e.g., 3-month-mean; Fig. 4b) and subseasonal (e.g., monthly; Supplementary Fig. 6) variations of SST and SLP anomalies associated with KE SSH. The close quantification of impacts of background state and ocean dynamics on air-sea turbulent heat exchanges over the KOE on subseasonal time scales might need to be revisited. Collectively, our results suggest that the summer KE SSH anomalies could be a critical component of predictability for the North Pacific seasonal air-sea coupled processes of earlier winter. The quantified oceanic and atmospheric components of earlier winter KOE surface heat fluxes show that the ocean-driven SST anomalies play a prime role in determining the ocean-atmosphere temperature difference, strengthening the anomalous heat transfer from the ocean to the atmosphere.
One of the merits of using the surface heat fluxes is that they can hint at slowly-evolving subsurface variations that are currently constrained by satellite observations (e.g., SSH) or Argo (e.g., ocean temperature). Given that the spatial patterns of KOE LHF/ SHF associated with KE SSH are characterized by quasi-stationary meanders related to the geographical location, we use projections of those patterns and extend their time series back to 1959 (see method for more details). The stationary spatial patterns of the KOE heat flux variability over >60 years in the observational records (i.e., 1959-2022) are confirmed in empirical orthogonal function (EOF) analysis of KOE LHF anomalies ( Supplementary Fig.  7), where the constructed KOE heat flux indices in this study (Fig. 7a) highly correspond to the most dominant variability (e.g., the leading EOF mode). Two-year running mean KOE LHF/SHF indices for the extended period in Fig. 7a exhibit that the covariation of LHF and SHF and their low-frequency components have noticeably increased since 1990. We attribute the increasing correlation of those two-year running mean KOE LHF and KOE SHF from 0.31 (before 1990) to 0.86 (after 1990) to their enhanced lowfrequency covariability because the correlation of monthly indices between KOE LHF and KOE SHF has barely changed from before (R = 0.72) to after (R = 0.75) 1990. We note that both the low-pass filtered KOE LHF and SHF variations have a prominent quasidecadal spectral peak (e.g., 8-12 years; Fig. 7b) that is identical to the dominant time scale of the KE SSH 2,21,22,56 . In particular, we detect a significant increase in the decadal variance (e.g.,~10year-period) of KOE heat fluxes since the mid-1980 in wavelet analysis of 2-year-running-mean KOE heat flux indices (Fig. 7c, d), which is confirmed both in the constructed indices of KOE LHF/ SHF and the leading principal components of KOE heat fluxes.
Consistent with our results that the decadal KOE heat flux variations are mainly driven by earlier winter air-sea heat exchange, the enhanced decadal fluctuations seem to be due to increasing variance of KOE LHF/SHF in OND and NDJ (Fig. 7e, f). It is important to note that increases in the seasonal variance of KOE LHF/SHF are detected only in the early cold season (i.e., OND and NDJ, Supplementary Fig. 8). We finally suggest that the increase in KOE heat flux variability in earlier winter can be mainly attributed to the enhanced oceanic variability (e.g., KOE SST) rather than atmospheric variance (e.g., wind speed) in Fig. 7g. All those increases in KOE variances exceeding~35% shown in Fig. 7e and g are highly significant (p < 0.01) based on a statistical significance test using the Monte Carlo approach (see method for details). The less significant increasing variance (Fig. 7f) of the leading time series of KOE LHF/SHF EOF might be due to strong stochastic components included in the dominant mode of surface heat fluxes. The observed changes in KOE variations in this study are consistent with recent observational and modeling studies that report significant changes in variance and time scale of the KE SSH, such as the stronger decadal KE variability during the recent period and under simulated future climate change 27,[57][58][59] . Those observed and projected enhanced decadal ocean variations in a warming climate suggest that the changes in earlier winter heat transport over the KOE might be a response to increasing anthropogenic forcing. However, without climate model experiments/projections and further analysis, we do not conclude that the increasing variance of KOE heat fluxes is solely due to external forcing because of the possibility of much lower frequency modulations (e.g., centennial time scales) of oceanic variations. Our study suggests that the increasing role of both physical (e.g., SSH) and thermodynamical (e.g., LHF/SHF) variations over the KOE region can be a key component of understanding the Pacific ocean-atmosphere coupled system as the conveyor of decadal ocean variability to air-sea coupling.
Given the increasing ocean-atmosphere coupled processes over the Northwestern Pacific on both seasonal and decadal timescales, one important and challenging question is how we can apply our findings and understandings to current climate predictions and climate change. Could the North Pacific WBC variations be important sources of seasonal-to-decadal climate predictability and change, especially for the Asia-Pacific-North America sector and variables? Further investigation of variations and changes in WBC thermal impact on midlatitude weather and climate, including hydrological cycles and extreme events, will be essential for real-time predictions and support for climate adaptation and mitigation.

Observational datasets
We use satellite SSH from the Copernicus Marine and Environment Monitoring Service (CMEMS), which has merged along-track satellite altimetry measurements since 1993 with a 1/12 o resolution. We use subsurface ocean temperatures from a reanalysis obtained from EN4 quality-controlled dataset 60 provided by Met Office Hadley Centre. To examine SST and atmospheric fields associated with KOE ocean variations, we use surface latent and sensible heat fluxes, sea level pressure, surface air temperature, 10 m zonal and meridional winds, and 10 m wind speed, which are taken from the European Centre for Medium-Range Weather Forecasts Reanalysis-5 (ERA5) 61 Fig. 9). All fields are obtained as monthly datasets (except for daily mooring data), and we use both monthly and seasonal mean (i.e., 3-month mean) to show decadal and seasonal evolutions of KOE variations, respectively. Monthly anomalies are computed by subtracting the climatological monthly means and the long-term linear trend using a least square method.

Removing tropical variability
Given the strong influence of tropical Pacific variability (e.g., ENSO) on midlatitude air-sea interaction, many previous studies of the North Pacific have attempted to remove tropical SST signals from the time series at all midlatitude grid points 2,62-64 . Following their approach, we apply and remove linear regressions onto the concurrent and lagged (1-3 month) ENSO signals from all the atmospheric and SST anomaly fields and WBC indices (e.g., KE SSH and KOE LHF/SHF time series) prior to analysis. The ENSO signal here is defined as the first three principal components of SST anomalies in the tropical Pacific between 12.5 o S and 12.5 o N following Frankignoul and Sennéchael 64 ( Supplementary Fig. 10). For example, each monthly anomaly at given longitude and latitude point X t ð Þ is replaced by X t ð Þ À aE 1 t À τ ð ÞÀbE 2 t À τ ð ÞÀcE 3 t À τ ð Þ, where E 1 t ð Þ; E 2 t ð Þ,and E 3 t ð Þ are those first three principal components of the tropical Pacific SST anomalies, τ indicates 0-3 months considering delayed response of atmospheric ENSO teleconnections, and a, b, and c are seasonally varying regression coefficients determined by least square fit for each variable at the grid point. The first mode (60% of the total variance) represents the traditional ENSO pattern, which is the tropical SST pattern regressed on Nino3.4 SST. The second mode (13.2% of the variance) represents the dipole SST contrast (i.e., east-west), and the third (4.4% of the variance) exhibits a narrow SST signature along the equator.

Definition of KE SSH and KOE LHF and SHF indices
We adopt an index of KE SSH variability introduced by Qiu et al. This index effectively describes dynamic properties of the Pacific WBC system, such as the latitude of the oceanic front and eastward jet, intensity of the KE recirculation gyre, and ocean heat content 9,15,33,38 . Spatial patterns of the KOE surface latent heat flux (LHF) and sensible heat flux (SHF) associated with the KE SSH are identified by regressing 3-month lagged LHF and SHF onto the KE SSH index. The 3-month lag is chosen based on the peak regression amplitude of lagged LHF/SHF anomalies regressed onto the KE SSH index (Fig. 2) 9 ). To obtain the corresponding KOE LHF/SHF timeseries, we project the spatial patterns of the KOE LHF and SHF onto monthly LHF and SHF anomalies, respectively, following Zhao and Di Lorenzo (2020) 65 who identify patterns and indices of ENSO precursor (note that we avoid using box-averages of KOE heat fluxes in order to detect geographical features and exclude unrelated intrinsic air-sea coupled processes). The domain for calculating the projections is the KOE region over the Northwestern Pacific (

Adjusted degrees of freedom
To conduct the significance test for correlations of decadal time scale and low-pass filtered climate indices (e.g., KE SSH index or 2 year-running mean time series), where the autocorrelation time scales are longer than~1.5 years, we use adjusted degrees of freedom following Bretherton et al. (1999) 66 : where N is the number of time steps, and ρ xx ðjÞ and ρ yy ðjÞ are the autocorrelation coefficient of time series x and y at lag j: Statistical significance tests using adjusted degrees of freedom have been widely used for investigation and assessment of low-frequency climate variations 67,68 .

Monte Carlo significance test
In this study, we show noticeable long-term changes in the earlier winter KOE variability for the extended observational record (i.e., 1959-2021). To determine if those annual changes (e.g., increasing/decreasing linear slope) are statistically significant, we use the Monte Carlo significance test with 30,000 random time series that contain the same number of time steps, mean, standard deviation, and autocorrelation as the observed KOE indices, via a first-order autoregressive model with no long-term trend. Using the obtained p-value based on distributions of slope changes of random sample variances, we determine the significance level that can reject a null hypothesis, which assumes that there are no changes in the linear slope (i.e., no long-term trends). For each synthetic time series, the running standard deviation is computed, along with the linear trend of that running standard deviation. The Monte Carlo sample percentiles of these trends provide tail probabilities for the null hypothesis of no forced trend, against which the observed trends can be compared.