Decadal variability of the Kuroshio Extension: the response of the jet to increased atmospheric resolution in a coupled ocean–atmosphere model

The Kuroshio Extension (KE) shifts between elongated and convoluted states on interannual to decadal time scales. The nature of this low frequency variability (LFV) is still under debate since it is known to be driven by intrinsic oceanic mechanisms, but it is also synchronized with the Pacific Decadal Oscillation (PDO). In this analysis we present the results from two present-climate coupled simulations performed with the CMCC-CM2 model under the CMIP6 HighResMIP protocol and differing only by their atmospheric component resolution. The impact of increased atmospheric resolution on the KE LFV is assessed inspecting several aspects: the KE bimodality, the large-scale variability and the air–sea interactions. The KE LFV and the teleconnection mechanism that connects the KE and the PDO are well captured by both configurations. However, higher atmospheric resolution favors the occurrence of the elongated state and leads to a more realistic PDO representation. Moreover, both simulations qualitatively capture the signatures of atmosphere-driven and ocean-driven regimes over the North Pacific Ocean, even if the higher resolution induces an excessively strong ocean–atmosphere coupling that leads to an overestimation of the air–sea feedbacks. This work highlights that the small scale atmospheric variability (resolution lower than 1°) does not substantially contribute to improve the realism of the KE LFV, but causes significant differences in the air–sea interaction over the KE region likely related to the strengthening of the coupling. The eddy-permitting ocean resolution shared by both configurations is likely responsible for the degree of realism exhibited by the simulated KE LFV in the two analyzed simulations.


Introduction
The KE is the eastward meandering jet formed by the convergence of the Kuroshio and Oyashio WBCs in the Northern Pacific Ocean and is known for exhibiting LFV over interannual and decadal time scales. The jet shifts between two main states: a quasi-stationary state characterized by two very intense anticyclonic meanders east of Japan with a strong zonal penetration (the so-called elongated state), and a meandering state characterized by a reduced zonal penetration in which the two anticyclonic meanders are absent (the so-called convoluted state). Two indices have been extensively used to highlight the KE bimodality (e.g., Chao 1984;Chen 2005, 2010;Pierini 2006Pierini , 2014Sasaki and Schneider 2011): the jet's latitudinal position <φ> and the path length LKE. Some degree of interdependency between these indices has been found. In particular, it has been observed that during the elongated state, a lower LKE is associated with a northward shifted jet, while during the convoluted state the jet shifts southward and the LKE substantially increases due to the formation of eddies following the dissipation of the two anticyclonic meanders east of Japan, in an inverse energy cascade process (e.g., Chen 2005, 2010;Gentile et al. 2018;Wyrtki et al. 1976;Yasuda et al. 1985;Yoon and Yasuda 1987).
Regarding the mechanisms underlying the KE LFV, two schools of thought with conflicting views exist. The debate revolves around the relative roles of intrinsic variability and external forcing in driving the quasi-decadal variability of the jet. The first view supports the idea that the KE LFV is an intrinsic mechanism led by the nonlinear interaction of recirculation, potential vorticity advection and eddies (e.g. Dijkstra and Ghil 2005;Ghil 2019;Jiang et al. 1995;Pierini 2006). Several studies (Dijkstra and Ghil 2005;Dijkstra 2005;Pierini 2006;Schmeits and Dijkstra 2001) show that the two regimes, in agreement with altimetry data (Qiu 2003;Chen 2005, 2010), can be reproduced using simple double-gyre idealized models forced with steadystate winds. These studies interpret the mechanism which leads the shift from the elongated state to the convoluted one as linked to the eastward propagation of potential vorticity anomalies from the western boundary current region (Cessi et al. 1987(Cessi et al. ,1990. In contrast, the second view supports the idea of a KE LFV induced by broad-scale atmospheric patterns of variability such as the Pacific Decadal Oscillation-PDO (e.g., Deser et al. 1999;Hogg et al. 2005;Mantua et al. 1997;Mantua 2002;Miller et al. 1994Miller et al. , 1998Nathan and Steven 2002;Qiu 2003;Qiu and Chen 2005). Taguchi et al. (2007), with the use of an eddy-resolving ocean general circulation model (OGCM), showed that the KE LFV is coherent with the large-scale atmospheric variability via the propagation of Rossby waves (RWs) from the East-Central Norther Pacific. In this context, it has also been shown that the frontal-scale variability contributes to changes in the KE jet speed which can be generated by intrinsic oceanic mechanisms, highlighting that while the linear Rossby wave theory explains the timing of the bimodality, nonlinear ocean dynamics plays an essential role in organizing the spatial frontal structure (Nonaka et al. 2006(Nonaka et al. , 2012Taguchi et al. 2007;Tsujino et al. 2006Tsujino et al. , 2013. Further studies, held by Pierini (2014) showed that the KE decadal bimodal cycle can be interpreted as a case of intrinsic variability paced by an external forcing in an ocean circulation model of the North Pacific under idealized forcings.
In order to better clarify the relative roles of internal and externally induced variability on the KE LFV, it is important to rely on realistic models of the ocean-atmosphere coupled system. Ma et al. (2015Ma et al. ( , 2016a analyzed the effect of air-sea coupling on the KE by comparing atmosphere-ocean coupled models at different resolutions. These results point to the relevance of explicitly resolving oceanic mesoscale eddies in order to faithfully reproduce the dynamics of, and the factors controlling, the KE jet. An increase of the spatial resolution in the oceanic component is relevant for the representation of ocean mesoscale turbulence and relative interaction with the atmosphere. On the other hand, the impact of the horizontal resolution of the atmospheric component on the main features of the KE jet and on its variability on quasi-decadal timescales has been relatively less inspected (e.g. Sakamoto et al. 2012;Delworth et al. 2012;Small et al. 2014).
In this context, by using coupled model configurations at different resolutions from the CMIP3 and CMIP5 archive (see Sect. 2.1), many authors inspected the role of ocean and atmosphere model resolutions on several aspects of the large-scale circulation, small-scale phenomena, extreme events, upwelling, oceanic eddies, sea-surface temperatures fronts, boundary currents and coastal currents (e.g. Sakamoto et al. 2012;Delworth et al. 2012;Small et al. 2014;Haarsma et al. 2016). Sakamoto et al. (2012), in particular, noted significant improvements in the degree of realism of the KE jet when the horizontal atmospheric resolution is increased in two coupled configurations: MIROC3m and MIROC3h (2.8125° and 1.125°, respectively), leaving the oceanic resolution unchanged (1.40625° × 0.56°-1.4°; zonal and meridional resolutions). The authors found that a finer horizontal resolution in the atmosphere introduces improvements in the orographic winds and on their effects on the ocean. In general terms, the benefits of increasing resolution on the degree of realism of the simulated phenomena is far from intuitive. In this context, Vanniere et al. (2019) showed that one cannot necessarily expect that all features of the simulated climate are in better agreement with the observations when run at higher resolution.
In this work we analyze the output of a fully coupled GCM, run under the HighResMIP experimental protocol (Haarsma et al. 2016). The latter has been specifically designed to inspect the effects of resolution on several key processes of the climate system. The specific goal of our analysis is to investigate the impact of model resolution on the representation of KE LFV, analyzed in terms of both frontal-scale and broad-scale variability, making a step forward into assessing the role of the horizontal atmospheric resolution on the degree of realism of the KE jet. The added value of the specific CMCC-CM2 implementation of the HighResMIP protocol does not only lie in the simple increment of the atmospheric resolution, but also in the matching between atmosphere and oceanic resolution, which is achieved in the VHR model configuration (where both ocean and atmosphere grids converge to the same 25 km grid-size). From other ongoing studies in the HighResMIP community (see some preliminary results reported in Tsartsali et al. 2019), this one-to-one gridpoint matching impacts on the strength of the ocean-atmosphere coupled interaction.
Both intrinsic and atmospheric processes play a fundamental role in the KE LFV, therefore it is necessary to inspect the role of the horizontal atmospheric model resolution in the representation of coupled ocean-atmosphere processes over this region. Atmospheric patterns of variability that interfer with the KE system, such as the PDO, are indeed likely affected by the atmospheric model resolution.
Being able to evaluate the benefits of advanced and wellevaluated high-resolution global climate models compared to their low resolution counterparts on the KE dynamical system, is fundamental for tuning the next generation models and define the best configurations able to simulate and predict regional climate with high reliability.
This work is organized as follows. In Sect. 2, the model, the observational datasets and the methodology applied in this study are described. Section 3 describes the results of the analyses, in terms of resolution impact on KE mean state, frontal-scale variability, large-scale variability and air-sea interactions. Finally, Sect. 4 concludes the paper.

Models and experimental setup
The model configurations used in this work are the Centro Euro-Mediterraneo sui Cambiamenti Climatici (CMCC) coupled atmosphere-ocean general circulation model CMCC-CM2-HR4 (High Resolution, HR; CMCC 2018a) and CMCC-CM2-VHR4 (Very-High Resolution, VHR; CMCC 2018b), which are part of the wider CMCC-CM2 family of models, documented in Cherchi et al. (2018). The HR and VHR configurations have been implemented and developed in the framework of the H2020 PRIMAVERA project (https ://www.prima vera-h2020 .eu), following the High Resolution Model Intercomparison Project (High-ResMIP) protocol for CMIP6. The models are composed of version 4 of the Community Atmosphere Model (CAM; Neale et al. 2013), NEMO3.6 for the ocean (Madec 2016), CICE4.0 for sea ice (Hunke et al. 2015), and version 4.5 of the Community Land Model (CLM, Oleson et al. 2013). Both HR and VHR configurations share the same eddy-permitting ocean resolution (ORCA025, i.e. 0.25°). The standard resolution configuration (HR) has an atmospheric resolution of 1° (64 km at 50° N), while the enhanced resolution (VHR) configuration has an atmospheric resolution of 0.25° (18 km at 50° N) (see Table 1 for a summary of the models' characteristics). In this respect, several studies showed the capability of ocean-only models to capture the KE mean field and variability under eddy-permitting configuration (0.25°; Pierini 2006Pierini , 2014Pierini , 2015Pierini et al. 2009).
Two 100-year control simulations of the present climate are conducted using fixed 1950s forcing conditions (referred to as "control-1950" in HighResMIP; Haarsma et al. 2016). Following the HighResMIP common protocol, the spin-up is initialized from the ocean at rest with temperature and salinity values corresponding to 1950 from the EN4 dataset (Good et al. 2013). After the spin-up, the models are further integrated for an additional 100 years with fixed 1950 forcing conditions. The analyses presented in this study focus on the 100-year model output, following the spin-up.
These simulations are the present climate equivalent of the CMIP6 DECK preindustrial runs (Eyring et al. 2016). Since the external forcings do not include interannual variability (i.e., no trends associated with observed changes in GHG concentrations, aerosol loadings, ozone and land use are included), only the LFV of intrinsic origin of the coupled system is represented in these control simulations. This set of model simulations provides an optimal framework to explore the sensitivity of the KE dynamics and its low frequency variability to changes in the model resolution. In the following analyses daily and monthly mean outputs of sea-surface height (SSH), ocean current velocities, sea-surface temperatures (SST) and turbulent (sensible and latent) heat fluxes at the air-sea interface are used. 1.2 × 10 -4 Background vertical eddy diffusivity (m 2 s −1 ) 1.2 × 10 -5 Isopycnal tracer diffusivity (m 2 s −1 ) 300 Eddy-induced velocity coefficient (m 2 s −1 ) n/a Albedo of snow on sea-ice No change Reference: Cherchi et al. (2018) 1 3

Observational data
For SSH, daily data on a 0.25° grid from the AVISO dataset (Archiving, Validation, and Interpretation of Satellite Oceanographic data; see AVISO webpage, http:// www.aviso .altim etry.fr, and the User Handbook) for the 1993-2017 period have been exploited. Gridded multisatellite geostrophic velocity (GSV) data, derived from Maps of Absolute Dynamic Topography (MADT), also provided by AVISO, were used to diagnose the surface flow pattern and the eddy kinetic energy (EKE) field. The latter is computed as follows: where ⋅ is the scalar product and ′ is the vector for the eddy component of the GSV field, defined as departures from the mean annual cycle. The EKE index is then obtained by averaging EKE over the upstream KE region [141°-153° E; 32°-38° N].
For the Ekman pumping velocity we use the 1° global monthly derived wind stress provided by NOAA (https :// coast watch .pfeg.noaa.gov) as follows: where is the wind stress, f is the Coriolis parameter and o is the reference density (assuming o = 1025 kg m −3 ). For SSTs and Sea Level Pressure (SLP), 1° monthly mean fields from the HadISST dataset (Rayner et al. 2003(Rayner et al. , 2016 have been used. The SSTs are used to diagnose the Pacific Decadal Oscillation (PDO) index from 1993 to 2017.
Finally, the ¼° daily SSTs from the NOAA OISST product (Reynolds et al. 2007) for years 1985-2013, combined with 1°, daily latent and sensible heat fluxes from the Woods Hole Oceanographic Institution objectively analyzed air-sea fluxes (OAFlux; Yu et al. 2008) were used to compute covariance maps and correlations between SST and surface turbulent heat fluxes over the KE region.

Climatology
Using daily SSH fields, two indices are computed to characterize the KE state: the jet's latitudinal position and path length. As already discussed in the introductory section, these indices are able to capture the KE bimodality from interannual to decadal timescales (e.g. Chen 2005, 2010;Pierini 2006Pierini , 2014. In the following, the two indices are computed in the upstream KE region (corresponding to sector 141°-153° E, shown in Fig. 1), where the highest jet variability is found due to the eddy-mean flow interaction. Downstream, the jet is still highly varying, but it is more influenced by the largescale (i.e., non-eddy-related) variability. Weekly values of the LKE and <φ> indices have been diagnosed, using a reference SSH isoline to identify the jet's axis and corresponding to the maximum meridional SSH gradient in the upstream KE region (0.6 and 0.7 m for HR and VHR respectively, while a 0.9 m value is used for the observations). No detrending filter has been applied to the indices since it is found that the residual drift (after the spin-up) explains a negligible fraction of the total variability of the KE system as described by LKE and <φ> indices. This was quantified by computing the ratio between the standard deviations of the total and detrended signals in both HR and VHR simulations, yielding O (1) values (not shown).
Simulated SSH climatology maps ( Fig. 1) show a latitudinal bias in the jet positions, with both HR and VHR showing a southward displacement with respect to observations of ~ 1° (see also the right panels of Fig. 3). Despite the bias, we can assess that both models are able to qualitatively capture the mean spatial structure of the SSH field in the Kuroshio and KE regions, with a strong anticyclonic meander south of Japan, a cyclonic circulation at the Kuroshio/ Oyashio confluence forming the KE, and two anticyclonic gyres east of Japan. Although both model configurations are able to reproduce some of the jet's climatological features, there are some differences that deserve to be commented. The SSH is clearly underestimated in both models, more in HR than VHR. However, HR better reproduces the mean field, with two, well distinguishable anticyclonic meanders east of Japan. In contrast, in VHR the upstream KE jet is more zonally elongated/less meandering compared to HR and observations. If the SSH mean field provides information about the KE mean features, the EKE mean field computed from GSV data, provides, to a good approximation, a measure of the eddy-driven variability. In Fig. 2 the simulated and observed EKE is shown.
The magnitude and location of the observed EKE signature is well captured by HR and VHR, but both models display a substantial underestimation of its zonal extension. However, HR shows a slightly less biased representation of the EKE pattern in terms of latitudinal extension, consistent with the more pronounced meandering character of the jet in HR, compared to VHR, the latter showing a more latitudinally confined and zonally stretched jet ( Fig. 1). In the next section the impact of resolution on the KE LFV is analyzed in detail.

Frontal-scale variability of the jet
The frontal-scale variability of the jet is characterized through the afore-mentioned LKE and <φ> indices. These are computed based on daily SSH field, which is previously low-pass filtered with a boxcar averaging using a 7-day smoothing window, in order to filter out the day-to-day variability. Figure 3 compares simulated and observed time series of the two indices. Note that model years bear no relation with the chronology shown for the AVISO data, the latter covering the 1993-2017 period. In order to better emphasize the variability over interannual and longer timescales, yearly filtered time series (obtained via a 1-year moving window averaging) are also shown in Fig. 3.
A simple visual inspection of the LKE time series (Fig. 3) reveals that the simulated and observed indices share a fairly consistent kind of variability, with sharp weekly-scale fluctuations embedded in a lower frequency, multi-annual scale envelope, with the jet length varying between an approximately 1500 km lower bound and higher values exceeding 3500 km.
The occurrence of specific regimes in the KE are better visualized by looking at the estimated probability density function (PDF) of the two indices. The well documented KE bimodality (e.g. Chen 2005, 2010) emerges from the AVISO data set (Figs. 4,5), displaying a more frequent, short LKE regime centered around a ~ 1500 km length, and a longer (but less frequent) LKE regime, with a ~ 2200 km length.
In order to allow a fair comparison between observations and models, accounting for the different lengths in the corresponding time series, we split the simulated LKE time series into 5 non-overlapping 20-year chunks and compare the simulated PDF, computed over each time segment, with the corresponding observational estimate. This comparison highlights the strong non-stationarity characterizing the modeled jet evolution, with different timeframes exhibiting substantial differences in their PDF distributions. The largest inter-decadal changes are featured by HR, alternating epochs with bimodal (1-20 and 80-100 time-frames in Fig. 4) and trimodal (40-60 and 60-80 timeframes) LKE distributions. VHR, on the other hand, reveals a closer adherence to the observed jet bimodality across most of its temporal evolution. However, the typical LKE values characterizing the model regimes and their relative frequencies undergo considerable changes across the different time segments. This, in turn, determines a substantial blurring of the pluri-modal jet behavior when the whole 100-years LKE time series is considered for both HR and VHR models, with the resulting PDFs showing a more unimodal distribution (Fig. 4, topand bottom-left panels). Both of them are skewed towards higher values, indicating that the occurrence of the elongated state (short LKE) is favoured over the convoluted state (long LKE), which can spread over a wider range of lengths (because of a high eddy activity during this state). This is in agreement with the observed distribution (in grey).
Regarding the <φ> index (Figs. 3, 5) both models show a clear (~ 1°) southward bias in the jet latitude. After splitting The low-pass filtered (1-year moving average) time series are also shown in red (green) for the LKE (mean latitude) index in the left (right) panels the 100-years time series into 20-year long, non-overlapping segments, a pluri-modal behavior emerges in both models, including sporadic hints of bimodality (notably, years 20-40 in HR and 60-80 in VHR, in Fig. 5). When considering the whole 100-year time series, the overall distribution appears to be largely Gaussian for HR, while some hints of bimodality are found in VHR.
To robustly assess the non-stationarity of the KE bimodality, the Kolmogorov-Smirnof test has been performed and applied (Kolmogorov 1933;Smirnof 1948) to the kernel density estimation of each 20-year PDF (Fig. 6).The null hypothesis (that states that the population is normally distributed) is rejected for every distribution and index with a 95% level of statistical significance. To confirm that this result is not timeframe-dependent we applied the same analysis to 24-year moving windows (consistent with the length of AVISO time series): again, each of the 76 distributions pass the statistical test.
The highly non-stationary nature of the jet is further highlighted by the continuous wavelet transform (Grossman et al. 1990) applied to the indices mentioned above (Fig. 7), thus providing a time-evolving picture of the simulated spectral features in the KE. To evaluate the significance of the wavelet transform we take advantage of the Matlab ASToolbox provided by Conraria and Soares (2011) and available at http://sites .googl e.com/site/aguia rconr aria/joana soare s-wavel ets. This test of significance is based on Monte Carlo simulations (Berkowitz and Kilian 2000) and confirms that the variability shown in Fig. 7 is statistically significant.
The wavelets highlight that the simulated KE properties, during the analyzed 100-year period, are highly nonstationary and that the mechanisms driving the variability act on a wide frequency spectrum. Periods of enhanced interannual-to-decadal variability in the evolution of the two indices emerge episodically. This non-stationarity of the system is higher in VHR than HR, with both configurations exhibiting the higher peaks of spectral energy around time scales of the order of 3-year and longer: these could be connected with a quasi-periodic atmospheric forcing.
To better illustrate the occurrence of bimodal behavior in the models' representation of the KE, we show in Figs. 8 and 9 the simulated weekly KE jet paths for a 20-year sample (corresponding to years 40-60, each panel corresponding to an individual model year) for HR and VHR, respectively. This time frame has been selected for the well distinguishable bimodality exhibited by the PDFs of the LKE and <φ> indices (Figs. 4,5), consistent with the concomitant emergence of interannual-to-decadal spectral features (Fig. 7). The sequence of yearly frames shows a consistent alternance of a convoluted (i.e., highly meandering, corresponding to high LKE values) and a zonally elongated (low LKE) state of the jet.
Several studies (e.g. Qiu and Chen 2010;Gentile et al. 2018) show that the LKE and <φ> indices are linked to each other during specific time intervals; in the elongated state (lower LKE) the jet is displaced northward, while during the convoluted state (higher LKE) the jet shifts southward. This led us to diagnose the Joint Probability Distribution (JPD) built upon the two analyzed jet indices, in order to identify the dominant regimes associated with the jet state in the bivariate <φ> -LKE phase space, for observations and models (Fig. 10). The KE jet bimodality emerges in the AVISObased JPD with two well distinguishable peaks: a higher frequency one corresponding to a northward jet with lower LKE (identifying the elongated state), and a less frequent  (Figs. 4, 5). The dotted black line indicates the kernel den-sity function for the simulated indices applied to the full timeseries (0-100 years) in HR and VHR. The thick black line indicates the AVISO kernel density function for each index one, southward-shifted and with a higher LKE, indicating a convoluted state (Fig. 10a).
The JPDs associated with the two model simulations span different areas of the <φ> -LKE space (Fig. 10b, c). HR covers a range of LKE values that is reasonably consistent with the observational estimates, with two local maxima (centered around 1700 and 2200 km) roughly corresponding to the ones detected in AVISO, but suffers from the previously discussed southward shift in the <φ> dimension. VHR, instead, under-represents the convoluted jet (high LKE) regime and further exacerbates (compared to HR) the southward displacement of the mean jet latitude.
Several studies highligth the relationship between the LKE and the <φ> index (e.g. Chen 2005, 2010;Pierini 2006Pierini , 2014, which is captured by AVISO in Fig. 10a: two main clusters representing the convoluted (higher length and lower latitude) and elongated (lower length and higher latitude) states are found. The significant correlation found in the observations between these two indices is about −29%. In contrast, in HR and VHR, significant correlations for each of the 76 distributions (24-year long) oscillate between −18 ÷ 33%, −28 ÷ 45%, respectively. External forcing and time-length are different in AVISO and in model simulations, therefore affecting their comparison. In particular, Fig. 10b, c show that the KE in the coupled simulations explore a bigger area of the latitude-LKE domain, which does not exclude the capability of the models to capture the elongated and the convoluted state but in different ranges of LKE and <φ> . It is interesting to note that, despite the detected biases, both models maintain a certain degree of jet bimodality, consistent with the observational estimates. However, these results also suggest that increasing the horizontal resolution of the atmospheric model alone does not contribute to mitigate the biases deriving from the use of a Fig. 7 Continuous wavelet transforms applied to the low-pass filtered indices: <φ> (left) and LKE (right) for HR (top) and VHR (bottom). The black contour designates the 5% significance level based on Monte Carlo simulations (Berkowitz and Kilian 2000). The cone of influence, which indicates the region affected by edge effects, is shown with a thick black line. The color code for power ranges from blue (low power) to red (high power). The white lines show the maxima of the undulations of the wavelet power spectrum 1 3 standard resolution model. A major impact of resolution is to favor the occurrence of zonally elongated (low LKE) jet regimes, compared to the convoluted jet regime. The models' ability in reproducing the observed PDO pattern is assessed in Fig. 11. Here, the PDO pattern is diagnosed through the leading EOF of SST anomalies in the [20° N-70° N] North Pacific domain, following Newman et al. (2016). Both simulated PDO patterns appear to be fairly consistent with observations. The only significant deviation is displayed by HR, featuring a spurious (negative) center of action off the eastern Japanese coast, which is absent in VHR. The regressed Sea Level Pressure anomaly (SLPa) on the PDO index for each dataset is also shown in Fig. 11. The differences between the regressed SLPa in the models and the observations can be expected since they refer to different forcing conditions and periods, while differences in the regressed SLPa in HR and VHR are the effect of the increased resolution of the atmospheric model on the PDO index. SLPa are very different between HR and VHR. In HR the regressed SLPa field displays a broad negative peak over the entire North Pacific and generally over the cold SSTa with a moderate amplitude (∼ −0.5 hPa) , while in VHR a well-defined minimum over the Gulf of Alaska with much larger amplitude (∼ −0.9 hPa) is found.

The influence of large-scale ocean-atmosphere variability on the jet
In order to understand how this distinct difference in the PDO-related SLP anomalies translates into the oceanic broad-scale features, it can be helpful to analyze the Ekman pumping velocity anomaly pattern, regressed on the PDO in every dataset. Figure 12 shows similar anomaly patterns south of the KE, with negative values that spread eastward from the east-central Pacific. These anomalies are stronger along the North American coasts in the observations, while in HR and VHR they reach the highest values in the east-central Pacific. North of the jet, every field has high positive anomalies along North America, that are more meridionally confined along the coasts in the observations, in contrast with HR and VHR where they extend seaward. Negative anomalies connect the central sector of the Pacific with the Japanese coasts in The PDO pattern discrepancies found when comparing VHR against HR are within the broad intra-model uncertainty assessed in CMIP models (see Fig. 1 in Wang and Miao 2018;and Fig. 8 in Newman et al. 2016). According to the comprehensive review presented in Newman et al. (2016) the PDO simulation uncertainty appears to be caused more by differences between models than by differences between realizations. Establishing what causes these modelto-model discrepancies is still under debate in the climate modelling community.
The link between large-scale variability in the extratropical North Pacific, as represented by the PDO, and the KE jet (LKE index) is illustrated in Fig. 13, showing the coherence between the PDO phase, the westward propagation of RWs and KE jet length variability (LKE). In order to better emphasize the RW signal, we applied a low-pass filter with a 1-year cut-off period to the indices and to the SSH anomaly field in order to reduce the monthly scale "noise". Here, the observed PDO index, SSH anomalies (Hovmöller diagram) and LKE index for the 1993-2017 period are shown. This multi-panel plot is an updated version of Qiu and Chen (2010) and captures a well-established teleconnection between changes in the broad-scale PDO variability mode and the frontal-scale variability of the KE jet, via the westward propagation of baroclinic RWs.
According to this picture, a positive PDO phase, associated with lower (higher) than normal SSH anomalies in the central North Pacific, leads, after propagation across the Pacific ocean basin, to an increase (decrease) in the KE path length, corresponding to a convoluted (elongated) state of the KE. As previously discussed, these LKE fluctuations are found to be correlated with a meridional shift of the jet (Qiu 2000), in turn associated with changes in the gyre circulation structure: during the PDO + (PDO −) phase, the Aleutian Low is stronger (weaker) than its normal, the Subtropical Gyre shifts southward (northward) and the negative (positive) SSH anomalies propagate westward, via baroclinic RWs, ultimately affecting the KE variability (e.g. Qiu 2000; Chen 2005, 2010;Di Lorenzo et al. 2008;Pierini 2014). Through this mechanism, SSH anomalies generated by PDO fluctuations in the eastern central Pacific (180°-200° E) pace the decadal scale variability of the KE.
Based on this interpretative framework, next we assess whether the same mechanism holds for HR and VHR models as well. Figure 14 shows the same diagnostics presented in Fig. 13, but for the two analyzed models. A number of similarities with the observational counterpart emerge, with This amplification is observed in many observational and modelling studies (e.g. Chen 2005, 2010;Wang and Wu 2019;Taguchi et al. 2005). Many possible mechanisms are associated with the amplification of the SSH anomalies in the Upstream KE region. First, the ocean internal variability is strong near the coast, which may contaminate the forced response. Another possibility is the interaction between baroclinic and barotropic modes (Taguchi et al. 2005).
In this respect, this amplification can be explained by the combined effect of the propagation of westward RWs and the origin of eastward RWs in the KE region (e.g. d 'Orgeville and Peltier 2009;Taguchi and Schneider 2014). Taguchi and Schneider (2014) used a 150-year coupled GCM integration and interpreted the eastward RWs as the counteraction of the jet to the arrival of westward propagating RWs in the KE region. These anomalies affect the mean temperature and salinity gradients, generating spiciness anomalies which are advected eastward by the mean currents. In response to these mechanisms the surface buoyancy flux generates higher baroclinic, eastward-propagating RWs.
A mechanism similar to the one suggested by Taguchi and Schneider (2014) may be occurring in both our models and could be related to the discontinuous propagation of the westward RWs in the central Pacific (Figs. 13, 14). In order to investigate this process, we have introduced a new index, defined as the Central Pacific (CP) index, computed as the SSH anomaly averaged over the [180°-200° E; 32°-34° N] lon-lat interval, roughly encompassing the PDO center of action in the two models. Given that to a negative PDO index corresponds a negative SSH anomaly, the CP index is defined as follows: This index is introduced to highlight the link between the variability in the east-central pacific (PDO center of action and therefore region of origin of the westward RWs) and the KE frontal variability. Since the PDO is a broader-pattern of variability, that includes many other phenomena in addiction to the westward RWs propagation from the east-central Pacific, we find that the CP index is more adapt to capture the teleconnection mechanism mentioned above. In this respect, when the PDO index and the CP index are not in phase (Figs. 13, 14) the RWs are not able to propagate continuously westward, suggesting that other phenomena may prevail. During these time intervals the PDO-related variability is overshadowed by other (non PDO-related) processes, including the eastward propagation of RWs. The Pearson correlations between the PDO and the CP indices is significant at the 95% level in AVISO (0.75), and in VHR and HR (0.60, 0.41 respectively). Therefore, an impact of the resolution is found in the reduction of the coherence between the CP and PDO indices. This is in agreement with Fig. 11 that shows a westward shift of the PDO center of action.
In order to better characterize the impact of the resolution, a more quantitative measure of the PDO-LKE connection is now obtained by computing the lagged correlation between the two indices for the 40-60 years time interval (Fig. 15) where the CP and PDO indices are highly correlated (~ 0.63 and 0.69 in HR and VHR respectively; p value ≤ 0.05 ). The PDO leads (lags) at positive (negative) in Fig. 16. In order to compute correlation from 0 to ± 10 years, the LKE and the PDO indices have been monthly averaged and a 1-year moving window filter has been applied. To compute the correlations coefficient in Fig. 16, we follow the Pearson's correlation and test the significance referring to the Student's t distribution (Emery and Thompson 2001), using the effective degrees of freedom to account for the self-correlation of the time series (Bretherton et al. 1999).
The correlations between the PDO and the LKE in these timeframes (Fig. 16) show a periodic correlation between the two signals in HR and VHR and highlight the differences in the teleconnection mechanisms between the PDO and the frontal variability of the jet, induced by changes in the resolution of the atmospheric model. When the PDO increases (decreases) the LKE increases (decreases) with a delay given by the time needed by the RWs to travel from the central Pacific (PDO center of action) to the KE region. The higher correlations (~ 0.8) are reached around ± 5 years, highlighting the interannual-to-decadal nature of this oscillation. In this respect, the negative/positive correlations at negative/positive time lags can be associated with the eastward/westward RWs propagation (e.g. d'Orgeville and Peltier 2009; Taguchi and Schneider 2014). Therefore, changes in the LKE are led by the PDO through westward RWs, and the PDO can be affected by spiciness anomalies that are advected eastward by the mean currents (eastward RWs).

Ocean-atmosphere interactions
The analysis presented in the previous sections is here complemented with an inspection of the air-sea interactions over a broad area surrounding the KE region. This aspect is relevant for both the frontal and broad-scale variability of the KE and can help interpreting the results presented in Sects. 3.2 and 3.3.
For this analysis the methodology illustrated in Bishop et al. (2017) has been used. In their work, the authors diagnose the lead-lag covariances between the SST and the turbulent (latent and sensible) heat fluxes (SHF)at the air-sea interface, and between SST tendency (SSTt) and SHF in several key areas of the world ocean, encompassing the major WBC systems. These relationships reveal whether the ocean-atmosphere energy exchange is driven by the atmospheric synoptic-scale variability (atmospheric weather) or by the oceanic internal variability associated with mesoscale eddies (oceanic weather). According to Bishop et al. (2017) observational analysis, the atmosphere-driven regime dominates in the open ocean (e.g., subtropical gyre areas), in contrast with WBC regions which appear dominated by the ocean-driven regime.
Here we adopt the same approach to verify whether the expected SST (and SSTt)-SHF lead-lag relations are faithfully reproduced by the two HR and VHR models. Specifically related to our case study, the following questions are addressed: are HR and VHR able to capture the shift from the oceanic-regime to the atmospheric-regime moving from the KE region to the open ocean? What is the impact of the horizontal atmospheric resolution in the ocean-atmosphere interplay?
The monthly lead-lag correlations have been computed between the SST and SHF and the SSTt and SHF for both the coupled configurations and the observations (see Sect. 2.2 for details on the observational products used in this analysis). For the heat fluxes we use the positive upward convention. Every field has been monthly averaged and 1-year moving window filtered. The lagged covariances are calculated as follows: The lagged covariance maps between the SST (SSTt) and the SHF (Figs. 17,18) in the models are in good agreement with the observational estimates but with some common behavior and differences. The SST (SSTt)-SHF covariances highlight a bias in the coupling strength, exhibited by both models when evaluated against the observations. This bias is further exacerbated by the enhanced resolution, as it is evident when comparing VHR with HR. In this respect, we added to Figs. 17 and 18 information about the significance of the covariances led by differences in the SST and SHF fields in HR and VHR. We find that differences between these fields in the models lead to significant covariances between the SST-SHF (p value ≤ 0.05) along the zonal extent of the KE, while they are negligible to the north and south of the jet. In contrast, the covariances between the SST tendency and SHF are significant (p value ≤ 0.05) upstream except for the near coast covariances and to the north and south of the jet. As done for Fig. 16, to test the significance we refer to the Student's t distribution, assuming the effecting degrees of freedom. These results confirm that horizontal atmospheric resolution higher than 1° causes significant differences in the air-sea interaction over the KE region, which can be related to the strengthening of the coupling in VHR. Generally speaking, the impact of model resolution on the air-sea "coupling strength" (Chelton et al. 2004) over western boundary current regions is the subject of ongoing work within the HighResMIP community. An extensive analysis of the coupling strength, as diagnosed via the pressure adjustment mechanism and vertical mixing mechanism in a set of HighResMIP simulations (Tsartsali et al. 2019) shows a substantial dependence of these mechanisms upon the models' grid resolution, with higher resolutions generally leading to a higher coupling strength. According to these analyses, both oceanic and atmospheric resolutions play a role on the coupling strength, but the strongest coupling is achieved when the ocean and atmosphere resolutions match, since under these circumstances the atmosphere does most effectively "perceive" the underlying oceanic variability features. Based on these analyses, we argue that the stronger air-sea interaction found in VHR, as compared to HR, reflects the above mentioned coupling strength sensitivity to atmospheric grid resolution, also noting that in VHR there is an actual matching between oceanic and atmospheric resolutions (25 km).
These biases don't affect the capability of the models to capture the main mechanisms that lead the variability in the KE region and in the open ocean. In fact, moving downstream from the WBC region, the variability shifts from an oceanic-driven to an atmospheric-driven regime. This is more evident in Fig. 19, which shows the lead-lag correlations between the SST (SSTt)-SHF computed in the WBC region and in the open ocean, choosing two points within and outside of the Kuroshio Extension (black diamond and green dot respectively; Figs. 17, 18). The lead-lag correlations computed in HR and VHR are in good agreement with the observations (Fig. 19).
In the open ocean the SST and the SHF are roughly in quadrature with each other around the zero lag and positive SHF. This ocean cooling is associated with negative SST tendency, identifying the atmospheric-driven regime. In the WBC region, the SHF acts damping the upper-ocean heat content anomalies generated by interior ocean processes, with the flux directly proportional to the SST itself, identifying the oceanic-driven regime.
To summarize, the lagged covariance maps (Figs. 17, 18) clearly show the ability of both model configurations to qualitatively reproduce the signatures of atmospheredriven and ocean-driven regimes over the open ocean and the WBC regions, respectively. The regions where the higher covariance values are observed, are reasonably well reproduced by the models. The zonal extension of the oceanicdriven regime is well captured by both the models, while its meridional extension is overestimated in VHR. In conclusion, we argue that the increased atmospheric resolution has a negligible effect on the oceanic variability in the KE region (where the intrinsic variability is always dominant) and that HR better captures the covariance relations compared to VHR, where an overly strong coupling strength leads to an overestimation of the air-sea feedbacks (Figs. 17,18).
Future investigations would certainly benefit from coupled simulations at finer ocean resolution, which would likely bring more consistent eddy dynamics and more

Final remarks
We presented an analysis of the KE decadal variability in two 100-year simulations of the present climate performed at standard and enhanced resolutions under the CMIP6 High-ResMIP common forcing protocol, with a state-of-the-art global climate model.
Both configurations display a reasonable degree of realism in reproducing the observed KE LFV. In particular, both systems exhibit multiple regimes in the evolution of the KE jet, featuring the alternation of convoluted and elongated jet length states, combined with northward/southward shifts in the jet positions. In turn, this variability appears to be affected by a strong non-stationarity, suggesting that the observed behavior might be representative of a subset of a potentially wider spectrum of possible KE jet states.
The main findings of this work are summarized below: • The KE's climatology is well captured by both configurations but enhancing the atmospheric resolution the jet appears to be less meandering. This is likely a rectification effect induced on the mean state of the jet by its frontal-scale variability (see the next bullet). • The frontal-scale variability in HR and VHR is consistent with the observational estimates, but an increased atmospheric resolution favors the occurrence of the elongated state, compared to the convoluted state. • A multi-modal behavior of the KE jet emerges in both models, including the intermittent emergence of bimodality. HR covers a range of LKE values consistent with the observational estimates but suffers from a bias in the latitudinal position of the jet. VHR, instead, underrepresents the convoluted jet regime and further exacerbates (compared to HR) the southward shift of the jet's position. Despite these biases, HR and VHR maintain a certain degree of jet bimodality, consistent with the observations.
The PDO acts as pacemaker of the KE LFV, in agreement with previous studies (e.g. Chen 2005, 2010;Pierini 2014), generating westward RWs that are able to pace the KE LFV. In this respect, we find that the CP index is more Fig. 18 Lagged covariances (− 1, 0, + 1 lag from left to right) between the sea surface temperature tendency SSTt and the surface heat fluxes SHF in the observations (top), HR (center) and VHR (bottom). Black triangle: WBC region; green dot: open ocean region. The gray shading identifies the areas where the differences between the SST and SHF in HR and VHR don't lead to significant covariances (p value > 0.05) adapt to capture the teleconnection mechanism mentioned above, since the PDO index includes a wider spectrum of phenomena. Moreover, when the PDO index and the CP index are not in phase the RWs are not able to propagate continuously westward, suggesting that other mechanisms may prevail (Sect. 3.3).
Increasing the atmospheric resolution leads to a more realistic PDO representation and the oceanic teleconnection mechanism provided by the propagation of baroclinic RWs is well captured in both models.
• HR and VHR are both able to qualitatively capture the signatures of atmosphere-driven and ocean-driven regimes of air-sea interaction over the open ocean and the WBC regions, respectively. However, higher resolution induces an overly strong air-sea coupling, leading to exceedingly large SST (and SSTt)-SHF covariances over the KE region.
This work highlights that, overall, the resolution increase of the atmospheric model does not substantially contribute to improve the realism of the KE LFV. The relatively marginal role played by the atmospheric model resolution may be explained by the primarily intrinsic nature of the KE LFV. The eddy-permitting ocean resolution shared by both HR and VHR configurations is likely responsible for the degree of realism exhibited by the simulated KE LFV in the two analyzed simulations.
On the other hand, even a substantial increase in the atmospheric resolution does not significantly improve the representation of the dominant mode of large-scale decadal variability in the extra-tropical North Pacific (i.e., the PDO; Sect. 3.3) while the biases in the patterns of air-sea interaction (Sect. 3.4) appear to be exacerbated by the resolution enhancement. These counter-intuitive aspects of model resolution enhancement are not novel in the recent literature (e.g., Docquier et al. 2019).
To conclude, both low and high-resolution simulations reproduce the observed bimodality of the jet and the oceanic teleconnection mechanism linking the PDO and the KE variability, via the propagation of baroclinic RWs. This essential mechanism is sufficiently well constrained under eddy-permitting ocean conditions, while the influence of the atmospheric resolution appears to be elusive, at this stage. Extending the present analysis to a multi-model multi-resolution ensemble of global climate simulations, benefitting from the whole set of HighResMIP integrations is a promising avenue for future investigations.
Funding The funding was received by PRIMAVERA EU Horizon 2020 Project (Grant no. 641727).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.