Atmospheric-Driven and Intrinsic Interannual-to-Decadal Variability in the Kuroshio Extension Jet and Eddy Activities

To investigate the influences of oceanic intrinsic/internal variability and its interannual-to-decadal modulations on the Kuroshio Extension (KE) jet speed and associated eddy activity, a ten-member ensemble integration of an eddy-resolving ocean general circulation model forced by the 1965–2016 atmospheric reanalysis is conducted. We found a distinct time–scale dependence in the ratio of forced and intrinsic variability of the KE jet speed. On the decadal time scale, the ratio of the magnitude of intrinsic variability to that of the atmospheric–driven variability is 0.73, suggesting it is largely atmospheric–driven. In contrast, on the interannual time scales, the KE jet speed has a large ensemble spread, indicating that it is strongly affected by intrinsic variability and has substantial uncertainty. For eddy activity, the ratios of atmospheric–driven and intrinsic variability also depend on the region. In the downstream KE [32°–38°N, 153°–165°E], variability in the atmospheric–driven eddy activity dominates (1.36 times) over the intrinsic variability on the decadal time scale, and is positively correlated with the current speed. Consistent with the westward propagation of atmospheric–driven jet speed anomalies shown by the ensemble mean, the eddy activity in the downstream KE region is correlated with the current speed variability in the central North Pacific 4 years earlier. This linkage is robust even for each ensemble member with the significant lagged correlation found in seven out of ten ensemble members as well as the ensemble mean (r = 0.59), suggesting the possibility of prediction of the eddy activity. In contrast, the eddy activity in the upstream KE [32°–38°N, 141°–153°E] shows a very large intrinsic and limited atmospheric–driven variability with a ratio of the former to the latter of 2.73. These results suggest that the intrinsic variability needs to be considered in the interannual variability of strong ocean jet. The dependence of these findings to the model specificities needs to be further explored.


INTRODUCTION
The Kuroshio is the western boundary current of the subtropical gyre in the North Pacific. After departing from the Japan coast at Cape Inubo at the eastern edge of the Honshu Island, it continues as an eastward free jet in the Kuroshio Extension (KE) and is associated with the strongest mesoscale eddy activity in the North Pacific Ocean (e.g., Qiu and Chen, 2005).
The Kuroshio and KE system transports heat from the tropics to the mid-latitude, and release huge amounts of heat to the atmosphere (about 1.7 × 10 8 W in the annual mean net surface heat flux from the KE region [30 • -40 • N, 141 • -165 • E], compared to 7.9 × 10 8 W in the whole Northern Hemisphere oceans to the north of 30 • N based on the third version of the Japanese Ocean Flux Data Sets with the Use of Remote Sensing Observations (J-OFURO3; Tomita et al., 2018) to impact the atmosphere aloft. Through such processes, the strong warm currents and associated sea surface temperature (SST) frontal structure affect the development of cyclones (Kuwano-Yoshida and Minobe, 2017;Hirata et al., 2018) and the large-scale atmospheric circulations (Minobe et al., 2008;Frankignoul et al., 2011;O'Reilly and Czaja, 2015), inducing basin-scale air-sea interactions (Qiu et al., 2014) extending to the tropical Pacific (Joh and Di Lorenzo, 2019). In addition, the decadal variability in the strength of the currents can affect SST, mixed layer depth, and further marine ecosystem in the Kuroshio Extension region (e.g., Nishikawa et al., 2011).
Eddies associated with the KE jet (Itoh and Yasuda, 2010;Sasaki and Minobe, 2015) are an important agent for water mass exchanges across the KE jet and associated oceanic frontal zones, and influence distributions of temperature (Sugimoto and Hanawa, 2011), potential vorticity (Qiu and Chen, 2006;Oka et al., 2015), nutrients (e.g., Sasai et al., 2019), and oxygen (Oka et al., 2015(Oka et al., , 2019. The eddy activity modifies also the KE jet itself (Ma et al., 2016). Further, eddies in the region may affect the ocean-to-atmosphere feedbacks (Ma et al., 2015). Because of these possible influences on the atmosphere and ocean, including geochemical variables, it is important to improve our understandings of variability in the KE and associated eddies.
Continuous satellite altimeter observations from the early 1990s have revealed that the KE jet has interannual-to-decadal variability in its strength, latitude, and its stability, which are tightly linked with the eddy activity (Qiu and Chen, 2005). The altimeter observations also unveiled that the interannual-todecadal variability in the KE jet strength is associated with the sea surface height anomalies that are generated in the central/eastern part of the North Pacific by wind-stress variability and propagate westward as Rossby waves (Qiu and Chen, 2005), indicating that the variability is atmospherically driven. On the other hand, idealized ocean model experiments suggest that strong western boundary currents have variability under steady atmospheric forcing (e.g., Dijkstra and Ghil, 2005). Similarly, a substantial interannual-to-decadal variability of the KE jet has been found in the layer and ocean general circulation models (OGCMs) that realistically represent the ocean circulation (Pierini, 2006;Taguchi et al., 2007;Sérazin et al., 2015) when driven by seasonally-varying atmospheric field without the interannual variability. To investigate atmospherically-forced and intrinsic oceanic variability and their non-linear interaction, ensemble OGCM simulations have been performed in recent years. The OCCIPUT project (OceaniC Chaos-ImPacts, strUcture, predicTability, Penduff et al., 2014) conducted and analyzed an ensemble of 50 eddy-permitting (1/4 • horizontal resolution) ocean model integrations subject to realistic atmospheric forcing for 56 years. This large ensemble revealed statistical properties of the atmospherically-forced and intrinsic oceanic variability around the globe (Penduff et al., 2019;Close et al., 2020) and in the Atlantic Meridional Overturning Circulation (AMOC; Leroux et al., 2018). South of the subpolar gyre, the AMOC simulated in 1/4 and 1/12 • horizontal resolution OGCMs showed no significant differences in the variance of intrinsic variability, validating the use of eddy-permitting models to investigate AMOC's intrinsic variability (Gregorio et al., 2015). Detailed structures of intrinsic variability of AMOC were further investigated by an eddy-resolving ensemble of multidecadal integrations for the North Atlantic Ocean (Jamet et al., 2019). Prominent intrinsic interannual variability in the KE jet was found in a 3-member ensemble of eddyresolving (1/10 • horizontal resolution), 18-year integration of a realistic OGCM forced by interannually varying atmospheric forcing under slightly different conditions (Nonaka et al., 2016). The result indicates that part of the KE jet fluctuations are caused intrinsically by non-linear interaction of the strong current and eddies. However, the limited ensemble size and integration period of Nonaka et al. (2016) experiment precluded a separation of interannual and decadal variability due to intrinsic ocean processes and atmospheric forcing, and pacing of intrinsic variability in the KE jet by atmospheric variability (Taguchi et al., 2007;Pierini, 2014).
The same question arises with respect to the interannual and decadal variability of ocean eddy activity. Interannual-todecadal variability in the KE jet strength, stability, and associated eddy activity induced by Rossby wave propagation (Qiu and Chen, 2005) suggests a wind-driven component of interannual variability in eddy activity. Meanwhile, as eddies are basically formed by oceanic dynamical instability, their activities could be inherently uncertain. In addition, the aforementioned intrinsic interannual variability in the KE jet speed (Pierini, 2014;Nonaka et al., 2016) is likely to induce intrinsic variability in eddy activity. Indeed, recent studies have shown that smaller horizontal-scale oceanic variability like mesoscale eddy activity is more strongly affected by oceanic intrinsic variability (Penduff et al., 2011;Sérazin et al., 2015). To further understand the interannual-todecadal variability and predictability/uncertainty of mesoscale eddy activity, it is necessary to clarify the potential importance of intrinsic variability.
The purpose of the present study is to investigate the importance of oceanic intrinsic variability in the KE jet strength and associated eddy activity on interannual-to-decadal time scales. This requires available computational resources be used to secure eddy-resolving horizontal resolution, albeit with a limited ensemble size. We conduct a 10 member ensemble of 52-year integrations of a realistic, eddy-resolving OGCM. In section 2, we introduce the OGCM and ensemble experiment, and show the long-term mean and its spread among ensemble members in section 3. We investigate the KE jet variability and eddy activity in sections 4 and 5, respectively. Section 6 and 7 are for the discussions and summary.

OFES2
The second version of the OGCM for the Earth Simulator (OFES2) interannual integration has been conducted for the period from 1958 to 2017 (Sasaki et al., 2018). For OFES2, we have modified the first version of OFES (Masumoto et al., 2004;Sasaki et al., 2008) in vertical mixing and tidal mixing parameterization, and coupled the ocean model with a sea-ice model (Komori et al., 2005). The integration has been performed with a horizontal resolution of 1/10 • in the model's near-global domain, from 76 • S to 76 • N. The topography of OFES2 is obtained from ETOPO1 (Amante and Eakins, 2009) with a maximum depth of 7,500 m. The model's vertical coordinate has 105 levels with 44 levels in the upper 300 m. The first version of OFES was integrated from 1950 (Sasaki et al., 2008) after a 50-year spin-up with climatological atmospheric field (Masumoto et al., 2004). OFES2's integration branched off this integration with the temperature and salinity initial conditions taken from the simulated fields of OFES on January 1st, 1958.
To estimate wind-stress at the sea surface, wind speeds relative to ocean currents are used in OFES2. The wind fields are obtained from the Japanese 55-year atmospheric reanalysis based surface dataset for driving ocean-sea-ice models (JRA-55do; Tsujino et al., 2018), in which wind fields are adjusted to the satelliteobserved wind speeds. As the satellites observe sea surface wind speeds relative to the surface ocean currents, the OGCM double counts the effect of the surface currents (Zhai and Greatbatch, 2007). Furthermore, the OGCM is driven by the atmospheric reanalysis field that does not include responses to the modeled SST variability, and could underestimate the western boundary currents and associated eddy activities (Renault et al., 2019(Renault et al., , 2020. Sasaki et al. (2020) shows further details of the model setup and comparisons between observations and modeled fields in OFES2. The model represents well observed fields with some exceptions, for example, a northward turn of the extension of the Gulf Stream in the North Atlantic.

Ensemble Experiment
In addition to the original OFES2 eddy-resolving interannual integration mentioned above, we conducted a 10-member ensemble integration of OFES2 for the period from the beginning of 1965 to the end of 2016 (OFES2_ensemble). Ensemble members differ in their initial conditions but are driven by the identical JRA-55do atmospheric reanalysis fields that were also used in the original OFES2 integration. Initial conditions are obtained by sampling the original OFES2 integration every 2 days over an 18 days period between January 3rd and 21st, 1965. Those fields were applied on January 1st, 1965 and then integrated to December 31st, 2016. From the slight differences in the initial conditions, differences among the ensemble members increase with time and saturate in about 5 years (not shown). We analyze the model output from 1970 to 2016. This method to generate the ensemble simulations is basically similar to that used by Jamet et al. (2019). As the differences developed under identical atmospheric variabilities, we attribute differences among the realizations to oceanic intrinsic variability. Temporal variability of the ensemble mean is assumed to largely result from atmospheric forcing, although the ensemble mean can include some components that are not atmospheric driven due to limited number of ensemble members.

Definition of KE Jet Speed and Eddy Kinetic Energy
In the present study, we define the KE jet speed as the current velocity at 2.5-m depth along its axis based on monthly-mean data. At each zonal grid point within the latitudinal band of 30 • -40 • N, we define the axis of the KE jet as the meridional grid point with the maximum current velocity. The occasional impact of mesoscale features on the jet axis is minimized by the use of a 10-degree longitude average and 13-month running mean (mentioned below).
We estimate the eddy kinetic energy (EKE) as (u 2 +v 2 )/2, where u and v are temporally high-pass filtered zonal and meridional velocity at 2.5-m depth, respectively, and u and v are calculated by subtracting the 13-month running mean fields from their original monthly time series. In the present study, we further apply 13-month running mean to the EKE and other time series to focus on the variability on interannual and longer time scales. To avoid long-term trends and/or multi-decadal variability, which could be induced by oceanic intrinsic variability (Sérazin et al., 2016) and can affect the analysis of time series, we removed the linear trends before the following analysis except for the analysis of the long-term mean fields for Figure 2E and relating discussion.

Observational Data
For observed data, we used the surface geostrophic currents distributed by the Copernicus Marine Environment Monitoring Service for comparison with the simulated data. The product identifier of the data is SEALEVEL_GLO_PHY_L4_REP_OBSERVATIONS_008_047. These currents are available on a 0.25 • × 0.25 • grid, and we use monthly-mean values estimated from the daily-mean data in this study. With the observational data, we estimate the KE jet speed and EKE in the same way as for the simulated data. We also use the indices of the North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al., 2008) and the Pacific Decadal Oscillation (PDO; Mantua et al., 1997) obtained from the web sites of the Georgia Institute of Technology 1 and Joint Institute for the Study of the Atmosphere and Ocean, University of Washington 2 , respectively. Data between 1970 and 2015 are used to examine correlations of the EKE in the KE region with these indices. J-OFURO3 data are obtained from the Asia-Pacific Data-Research Center at the International Pacific Research Center, University of Hawaii 3 .

Estimation of Magnitudes of Atmospheric-Forced and Intrinsic Variabilities
In the following analyses, we hypothesize that the variability of the ensemble mean f of variable f represents the ocean response to time dependent atmospheric forcing, and that deviations f − f of each ensemble member from the ensemble mean are caused by intrinsic oceanic variability. Here, the angled brackets represent the ensemble mean, f = 1 M M m=1 f , and an overbar temporal average,f = 1 T T t=1 f , where M is the ensemble size and T is number of months for the analysis period. We define the ensemble spread as the standard deviation among ensemble members from the ensemble mean: Following Rowell et al. (1995), see also Leroux et al., 2018), the variance of intrinsic variability is estimated by: the variance of the ensemble mean is: and the variance of the response to variable atmospheric forcing is: In σ 2 atm , the term of − 1 M σ 2 int is added to take into account the influence of the limited number of ensemble members. The magnitudes of the intrinsic and atmospheric-driven variability are characterized by σ int , and σ atm .

Global Distributions
We first examine the influence of oceanic intrinsic variability on long-term, 30-year (1986-2015), mean fields. Figure   Frontiers in Marine Science | www.frontiersin.org indicates that even in the 30-year mean current speed, there are spreads among the ensemble members in the western boundary currents and the Antarctic Circumpolar Current regions. Corresponding ensemble spreads are also found in the same regions in the 30-year mean sea surface temperature (SST) field ( Figure 1B). While the modeled SST field tends to be restored to observed surface air temperature field through the estimation of surface heat fluxes based on the bulk formulae, SSTs are strongly influenced by oceanic dynamics in the strong current regions and can have the ensemble spreads associated with the spreads of the current field. The distribution of the spread in the long-term mean fields is consistent with regions of high intrinsic variability found in the previous studies based on climatological atmospheric forcing (Nonaka et al., 2016) and realistic atmospheric forcing (Close et al., 2020). Such intrinsic variability induces differences among the ensemble members even in the long-term mean fields. However, due to limited number of 10 ensemble members, estimates of the ensemble spread have uncertainties corresponding to the particular ensemble size as discussed below.

The Kuroshio Extension Current
The long-term (30-year) and ensemble mean of the surface current in the North Pacific (Figure 2) indicates that the mean current distributions are well represented in the model, although speeds are lower than the observations (Figure 2A), even if we compare the observed data to the ensemble member with the strongest current in the KE region ( Figure 2B) or change the averaging interval to the overlapping period from 1993 to 2016 (not show). The surface boundary condition for the wind stress based on the wind speeds relative to ocean currents could be a reason for the weaker surface current.
As discussed above, the ensemble spreads of the long-term mean surface current are large along the Kuroshio and KE, and extends to around 170 • E ( Figure 2D). While absolute values of the ensemble spreads are limited to less than 5 cm s −1 , except for the Kuroshio large meander region, they exceed the atmospherically-forced variability in different 30-year averages of the ensemble mean. Indeed, the ensemble spread ( Figure 2D) is generally larger than the differences between the two 30year averaged ensemble means for distinct averaging intervals such as  and   (Figure 2E), the most separated 30-year means in the period of 1970-2016. Mean differences among four 30-year averaged ensemble means for the averaging intervals starting in 1971, 1976, 1981, and 1986 are further small (not show), and thus smaller than the ensemble spread.
Oceanic internal dynamics leave a clear difference in the longterm mean jet structures averaged from 160 • E to 170 • E, an area of large ensemble spreads and weaker mean current speed of the KE. Figure 3 shows that the mean KE has substantially different structures among the ensemble members: some members have rather broad single jet, while other members have double or even triple jet structures (Figures 3A-C). Consistently, the relatively large (about 1 cm s −1 ) ensemble spreads extend vertically to about 300 m (Figure 3E), the typical depth of the jet structure (E) Absolute value of difference between two long-term mean simulated surface current fields for the periods of  and . Unit is cm s −1 .
found in the ensemble mean ( Figure 3D). Ensemble spreads of temperature are found in the thermocline layer below the enhanced spread of the current speed, consistent with the thermal wind relationship (Figure 3F). This suggests that the observed jet structures are difficult to reproduce in eddy-resolving OGCMs in this particular region due to its substantial uncertainty. Figure 4A depicts the time series of the KE jet speed for each ensemble member and ensemble mean for the integration period, excluding the first 5 years when the difference among the ensemble members develops. The ensemble spread has a substantial amplitude compared to the total variability of the ensemble average low pass filtered with a 13-month running mean. Indeed, the magnitude of its intrinsic variability, σ int (4.1 cm s −1 ), is comparable to the magnitude of atmospheric-driven variability, σ atm (4.4 cm s −1 ), and their ratio is 0.93. However, on decadal time scales, the spread among the ensemble members is rather limited, and all members have similar variations with peaks around 1974, 1980, 1990, and 2005. To focus on the decadal time scales, the KE jet speed time series is smoothed by a 37-month running mean: in this time series, which are dominated by decadal variability (not show), the ratio of magnitude of intrinsic variability σ int (2.8 cm s −1 ) to that of atmospheric-driven one σ atm (3.8 cm s −1 ), 0.73, is lower than the counterpart, 0.93, for the variability after 13-month running mean. In contrast, on interannual time scales, the spread among the ensemble members is substantial compared to the variability of the ensemble mean. For the band-pass filtered time series of [13-month running mean]−[37month running mean], which are dominated by interannual variability (not shown), the magnitude of intrinsic variability σ int (2.7 cm s −1 ) is larger than that of atmospheric-driven one σ atm (1.4 cm s −1 ), and the ratio of it is 1.92. In short, in the ensemble members, the KE jet speed exhibits on the interannual and longer time-scale a similar magnitude of atmosphericdriven and intrinsic variability, but its decadal and interannual components are dominated by atmospheric-driven and intrinsic variability, respectively.

Interannual-to-Decadal Variability in the KE Jet Speed
Comparing the time series of the ensemble mean and ensemble spread confirms above-mentioned features ( Figure 4B). Although the spread tends to be high (low) around decadal local maxima (minima) of the ensemble mean in 1974 and 1990 (1985 and 1997), such a relationship is not present in other peaks of the ensemble mean, and we do not find clear relationships between the variability in the ensemble mean and spread. It is indeed very clear that different time scales are at play: the ensemble mean is dominated by decadal time scale, and variability in the spread does not show clear decadal signal and is dominated by interannual time-scale. This is consistent with the rather limited ensemble spread on the decadal time scale and means that the decadal variability in the KE jet speed is more atmospheric-driven, and the KE jet speed has a large uncertainty on the interannual time scale. The thin red curve in Figure 4B indicates robustness of the variance among the ensemble members following the method of Leroux et al. (2018, Eq. A1). Due to the limited number of ensemble members in the present study, the uncertainty is very large compared to the ensemble spread itself. Therefore, we need to consider the ensemble spread as a result of the limited number of samples and to interpret the results with sufficient care. Comparison of the model results with observations (blue thick curve) consistent with the time-scale dependent uncertainty of the KE jet speed ( Figure 4A). As discussed with Figure 2, the time mean KE jet speed is lower in the model (93.7 cm s −1 in July 1993-June 2016, the common period for the model and observation) than the observation (110.4 cm s −1 ). Meanwhile, the ensemble mean represents very well the observed variability on the decadal time scale. Their correlation coefficient r = 0.76 (0.90) for the time series applied 13-month (37-month) running mean, 95 (99) % significant with the effective degree of freedom N = 7 (5) estimated following Metz (1991). In contrast, the band-pass filtered ensemble mean and observation, which are dominated by the interannual time scale, are unrelated (r = 0.17, not significant at 90% level with N = 23).

Eddy Activity in the KE Region
In addition to the simulated KE jet being slightly weaker than the observed (Figure 2), EKE in the model is underestimated compared to the observations ( Figure 5A) even in the ensemble member with the maximum EKE level in the high EKE region ( Figure 5B). However, the distribution of EKE is similar in the ensemble member with the minimum EKE in the region ( Figure 5C) and also in the ensemble mean (Figure 5D), and the model represents the long-term mean horizontal distributions well. Then, we investigate the temporal variability of simulated  Figure 5D).

KE Downstream
In contrast with the KE upstream region (Figure 6A) that will be discussed later, EKE time series in ensemble members in the downstream region show only small differences among the ensemble members especially on the decadal time scales (Figure 6B). Indeed, the atmospheric-driven variability in the eddy activity dominates (1.36 times) over the intrinsic variability on the decadal time scales (37-month running mean filtered time series), while the ratio is 0.34 on an interannual time scale (bandpass filtered time series). Also, the ensemble mean corresponds well with the observed time series except for the last several years (r = 0.73 for 1994-2010, 99% significant, N = 12, r = 0.41 for 1994-2015, not significant). The reasons for the difference between the observation and the simulation after 2010 are unclear at this stage. This means that in this region, at least in this particular model, the modulation of the eddy activity has a substantial atmospherically-driven component and thus, some potential predictability.
To investigate the mechanisms for the variability in eddy activities, we compare the longitude-time section of anomalies of the ensemble member 2 of EKE and current speed meridionally averaged from 32 • to 38 • N in Figure 7A. Both anomalies tend to propagate westward from the central North Pacific with an intensification of the EKE amplitude toward the western boundary region, and tend to co-vary in the KE downstream region. Indeed, the time series of area means in the KE downstream region (Figure 7B) are synchronized with each other: high current speed is accompanied by high EKE (r = 0.79, 98% significant, N = 8). For each ensemble member, the corresponding correlation varies from r = 0.79 to 0.93, but all of them are significant at 98% or higher level. The ensemble member 2 shown in Figures 7A,B is the member that has the lowest correlation between the time series shown in Figure 7B, among the all members. The relationship between the current speed and EKE is also suggested in the observations (Figure 7C, r = 0.89, 99% significant, N = 7). The westward propagation of the signals ( Figure 7D) and correlation between the area averaged EKE and current speed in the KE downstream region ( Figure 7E) are also found in the ensemble mean (r = 0.93, 99% significant, N = 6), indicating that the variabilities include the atmospheric-driven component.
As the current speed anomalies include the atmosphericdriven component and tend to propagate westward (Figure 7D), the propagation of anomalies could provide predictability for EKE in the KE downstream region. To explore this hypothesis, we first attribute changes of EKE to atmospheric forcing by examining the ensemble mean (Figures 8A-D). Predictability and comparison with the single observed realization further require estimation of the signal to noise ratio, estimated from the relative variances due to atmospheric forcing and intrinsic variations, and illustrated by the correlations for individual ensemble members (Figures 8E,F). In the lagged correlation maps, a high correlation region is found with current speeds even in 4 years earlier in the central North Pacific (Figure 8C), and then propagates westward for a shorter lead-time (Figures 8A,B). This feature is consistent with Taguchi et al. (2010) who showed that the low-frequency modulation of EKE in the KE region is associated with the incoming, westward-propagating Rossby wave signals. The time series of area averaged current speed in the central North Pacific (34 • -36 • N, 175 • W-165 • W) correlates well with the eddy activities in the KE downstream region in 4 years later (r = 0.59, 90% significant, N = 12), indicating a potential of predictability of the ensemble mean eddy activities in the region (Figure 8D). While Figure 8B suggests that the current in the same region also correlates with the eddy activities in the KE downstream 2 years later, the correlation is slightly lower (r = 0.50) and not significant. Also, the westward propagating signals are not found in the 6-year lead correlation map (not shown).
In the actual prediction, however, we need to consider the lagged correlation for each ensemble member rather than the ensemble mean as real observation is equivalent to one realization of the ensembles. For each ensemble member, the lagged correlation between the area averaged current speed in the central North Pacific and EKE in the KE downstream region 4 years later varies from r = 0.264 (ensemble member 6, Figure 8F) to 0.602 (ensemble member 3, Figure 8E) and their average is r = 0.446. While seven of ten members show the lagged correlation statistically significant at 90% level, the correlations are not significant for three members. (It should be noted that as the ensemble mean has smoother time series, the effective degree of freedom is smaller than that for each ensemble member, and correlation coefficients at 90% significance higher for the ensemble mean than for each ensemble member). This means that while the ensemble mean shows a potential of a 4-year lead prediction of eddy activity in the KE downstream region, it is not always expected for each ensemble member due to the intrinsic variability and, thus for real observation.
Indeed, a similar 4-year lagged correlation between the current speed in the central North Pacific (40 • -42 • N, 160 • W-150 • W) and EKE in the KE downstream region is also found in the observed data (r = 0.60), but it is not statistically significant (not shown). While this result can be due to intrinsic variability included in the observation, the less robust statistical relation could be also due to the limited length of observed time series. Therefore, to determine the potential predictability in EKE in the KE downstream region with several years lead-time in the real North Pacific Ocean has to await a longer observed record.

KE Upstream
In contrast to the downstream region, the time series in the KE upstream region of EKE for each ensemble member ( Figure 6A) shows substantial differences and there is almost no common temporal variability among the ensemble members. Indeed, the ratio of the magnitude of intrinsic variability of EKE to that of the atmospheric-driven variability is 2.73 in this region. Hence, at least in this particular model, the variability in the eddy activity in the region is strongly uncertain, and consistently, the time series based on observation does not agree with the model ensemble mean (r = −0.16, not significant). Nevertheless, the observed time series is in the range of the ensemble spread except for some years in the mid and late 1990s, implying that the observation is basically within the simulated uncertainty. We will discuss the high uncertainty of EKE in this region and the influence of model biases on the uncertainty in section "Discussion About the Simulated Eke in the KE Upstream Region."

Map of Potential Predictability
In the above analyses, we have investigated uncertainty/potential predictability of the eddy activity focusing on the two particular regions, the KE upstream and downstream regions defined by Qiu and Chen (2005). To examine horizontal distributions of this property in the whole KE region, we further plot the horizontal map of potential predictability for the interannual and longer time scale variability (Figure 9). Here, we define the potential predictability (PP) as follows: PP = σ 2 atm / σ 2 atm + σ 2 int (Rowell et al., 1995;Sugi et al., 1997). PP represents the ratio of deterministic, atmospheric-driven, variance to the total variance, considering the variability in the ensemble mean is deterministic.
In Figure 9, we plot the square root of PP estimated for 5 •longitude and 5 • -latitude mean eddy activities after applying a 13-month running mean to focus on the interannual and longer variability. The linear trends are also removed. The plot shows that PP is generally higher in the downstream side with peaks around 160 • E to the north of the KE jet, and has a local minimum around the KE jet (∼37 • N) in the further downstream region. It is clear that PP is low in the upstream of the KE jet around its axis, consistent with the difference found in the specific KE upstream and downstream regions in the previous subsections.  Comparison with the long-term mean horizontal distribution of eddy activity (Figure 5) implies that PP tends to be lower with higher eddy activity in the KE region.

DISCUSSION OF THE SIMULATED EKE IN THE KE UPSTREAM REGION
The relationship between EKE and the current speed found in the KE downstream region can be explained by dynamical stability. Given the geostrophic/thermal wind relationship with weak currents at depth, the stronger surface zonal current is associated with a stronger meridional temperature/density gradient. Stronger meridional temperature gradients are more baroclinically unstable and induce higher eddy activity. Stronger meridional shear associated with the stronger zonal currents favors barotropic instability.
The relationship we identified in the KE downstream region is, however, opposite to what Qiu and Chen (2005) found based on satellite observations in the KE upstream region. They show that a high (low) eddy activity is associated with a low (high) current speed on the decadal time scale. In the KE upstream region, as found in Figure 6A, the model does not capture the observed EKE variability. Furthermore, the observed relationship between the current speed and EKE variability is not clearly reproduced even in each ensemble member (only two of ten members show significant negative correlation between them consistent with the observation) as well as in the ensemble mean. So, the discrepancy in the region could be a model deficiency. Alternatively, this may imply that EKE activity in the KE upstream region is highly uncertain in the real ocean, so that the observed time series cannot be reproduced in ocean GCMs. Sugimoto and Hanawa (2012) show from observational data that when the Kuroshio takes its offshore non-large meander path and passes the southern part (32 • -33 • N) of the Izu Ridge (around 140 • E), eddy activity tends to be high in the upstream KE region. While OFES2 has a bias in the path of Kuroshio south of Japan, and tends to have the large meander path more often compared to the observations as suggested in Figure 2, we have confirmed that the model tends to take a near shore path at the Izu Ridge more often than in the observation (not shown). Based on the results of Sugimoto and Hanawa (2012), this means that the bias of the model Kuroshio path does not enhance the eddy activity in the upstream KE region. Also, as discussed with Figure 5, the model has slightly weaker EKE in the KE region compared to the observation, and the high uncertainty in the KE upstream region is not due to too strong eddy activity in the region.
The uncertainty in EKE in the upstream KE region can be further explored by correlating across all the ensemble members the intrinsic (unforced) component of SSH anomalies onto the intrinsic component of the EKE in the upstream KE region (Figure 10). Positive EKE deviations in the upstream KE region that exceed the ensemble mean are associated with a negative (positive) SSH deviations from the ensemble mean to the south (north) of the upstream KE jet ( Figure 10A). Hence, as the intrinsic variability, the stronger EKE is associated with a weaker southern recirculation gyre and a weaker SSH meridional gradient, while there is a possibility that the weaker recirculation and meridional SSH gradient results from the stronger EKE. This property is consistent with observation (Qiu and Chen, 2005), although the observed data includes both the atmosphericforced and intrinsic components. Interestingly, stronger EKE in the upstream KE region is also associated with positive SSH deviations in the Kuroshio south of Japan, and the lagged correlation maps (Figures 10B,C) show that the signal appears first in the region and then expands into the upstream KE region. The positive SSH deviations to the north of the mean Kuroshio meandered path (depicted by the contours) tends to straighten the Kuroshio path and this leads to stronger EKE in the upstream KE region 2 years later. The analysis suggests that the selfsustained intrinsic path variations of Kuroshio south of Japan (e.g., Qiu and Miao, 2000) may have a down-stream influence on the uncertainty of the KE's southern recirculation gyre and the EKE in the upstream KE region.

SUMMARY
To investigate the possible influence of oceanic intrinsic variability and its interannual-to-decadal modulations in the KE jet speed and associated eddy activities, we conducted a 10-member ensemble, 52-year integration of an eddy-resolving OGCM with time-varying surface forcing and slightly different initial conditions. We show that, on the decadal time scales, variability in the KE jet speed has a limited ensemble spread, suggesting that the decadal jet speed variability is mostly atmospheric-driven. In contrast, on interannual time scales, the KE jet speed has a large ensemble spread, indicating that it is strongly affected by intrinsic variability, has substantial uncertainty, and is difficult to predict. The present study confirms with a longer integration results of Nonaka et al. (2016). This time-scale dependence is similar to that found in the Atlantic Meridional Overturning Circulation (Jamet et al., 2019). While Pierini (2006) suggests the importance of intrinsic variability to decadal KE variability, the present study suggests that the intrinsic variability is more important on the interannual time scale. The pacing of intrinsic variability in the KE jet by atmospheric variability proposed by previous studies (Taguchi et al., 2007;Pierini, 2014) could not be identified in the relation between the intrinsic and atmospheric-forced variability of the upstream KE jet speed. Further exploration is necessary.
In the KE upstream region, observations show that eddy activity is high (low) when the KE jet is weak (strong) (Qiu and Chen, 2005). However, in our simulations, the eddy activity in the KE upstream region shows high uncertainty and limited atmospheric-driven variability. In contrast, the eddy activities in the KE downstream region have a small ensemble spread especially on the decadal time scale, indicating they are predominantly atmospheric-driven. Further, the eddy activity is highly correlated with the surface current speed averaged in the same region. This relationship is consistent with the enhanced baroclinic instability associated with stronger currents and meridional density gradients. Through westward propagation of surface current variability signals shown by the lagged correlation maps based on the ensemble mean field, which can represent atmospheric-driven potentially predictable component of variability, the eddy activity in the KE downstream region correlates well (r = 0.59, 90% significant) with the ensemble mean surface current speed variability in the central North Pacific 4 years earlier. Even for each ensemble member, the significant lagged correlation is also found in seven out of ten ensemble members, suggesting a possibility of prediction of the eddy activity in the region based on westward propagation of wind-driven current speed anomalies. Consistent with this, the ensemble mean of eddy activity in the KE downstream has lagged correlation with the index of NPGO, which leads the KE jet variability with several years lag. Their correlation is r = 0.41 (90% significant, N = 22) when the NPGO index leads a 39 months with applying 13-month running mean. Interestingly, the ensemble mean of eddy activity in the KE upstream correlates with PDO (r = 0.42, 99% significant, N = 48), when the PDO index leads 13 months. As we discussed above, however, the wind-driven component is limited in the KE upstream region and this correlation is not found for each ensemble member at least in our simulation due to the large uncertainty caused by intrinsic variability.
While it has been widely recognized in the atmospheric science that observations are just one realization of possible fields, the oceanography community is only beginning to explore this aspect (e.g., Penduff et al., 2018), except for the fields directly relating to eddies, and should consider oceanic fields in the similar way. There are a few emerging projects that have pointed out the importance of oceanic intrinsic variability (see Introduction). Among them, statistical properties of the atmospheric-forced and intrinsic oceanic variability have been revealed globally with a large ensemble OGCM simulations performed under the OCCIPUT project (e.g., Penduff et al., 2019;Close et al., 2020). In the present study, as we used more computational resources for high horizontal resolutions to represent the KE jet structure and to resolve mesoscale eddies, the number of ensemble members (10) is limited compared to the OCCIPUT experiment (50). While this choice precluded us from the robust estimates of ensemble spreads and Gaussianity of distribution (Penduff et al., 2019), our higher-resolution ensemble simulations complement the OCCIPUT simulations. By better resolving the eddy activity and jet structures, we have documented the atmospheric-forced and intrinsic variability on the interannual-to-decadal time scales in the KE region. The present results confirm that, even on the interannual time-scale that is longer than a typical eddy lifetime, it is necessary to consider uncertainty in variability of the strong jets. Nevertheless, we found possibility of predictability in the eddy activity in the KE downstream region. Since intrinsic variability is likely model dependent, these results should be explored with other models or using multi-model comparisons as well as larger number of ensemble members.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.