Imprint of chaotic ocean variability on transports in the southwestern Pacific at interannual timescales

. The southwestern Paciﬁc Ocean sits at a bifurcation where southern subtropical waters are redistributed equatorward and poleward by different ocean currents. The processes governing the interannual variability of these currents are not completely understood. This issue is investigated using a probabilistic modeling strategy that allows disentangling the atmospherically forced deterministic ocean variability and the chaotic intrinsic ocean variability. A large ensemble of 50 simulations performed with the same ocean general circulation model (OGCM) driven by the same realistic atmospheric forcing and only differing by a small initial perturbation is analyzed over 1980–2015. Our results show that, in the southwestern Paciﬁc, the interannual variability of the transports is strongly dominated by chaotic ocean variability south of 20 ◦ S. In the tropics, while the interannual variability of transports and eddy kinetic energy modulation are largely deterministic and explained by the El Niño–Southern Oscillation (ENSO), ocean nonlinear processes still explain 10 % to 20 % of their interannual variance at large scale. Regions of strong chaotic variance generally coincide with regions of high mesoscale activity, suggesting that a spontaneous inverse cascade is at work from the mesoscale toward lower frequencies and larger scales. The spatiotemporal features of the low-frequency oceanic chaotic variability are complex but spatially coherent within certain regions. In the Subtropical Countercurrent area, they appear as interannually varying, zonally elongated alternating current structures, while in the EAC (East Australian Current) region, they are eddy-shaped. Given this strong imprint of large-scale chaotic oceanic ﬂuctuations, our results question the attribution of interannual variability to the atmospheric forcing in the region from pointwise observations and one-member simulations.

Mesoscale eddies, whose evolution is largely chaotic, are ubiquitous in the region south of about 20°S (Keppler 80 et al. 2018), and the currents variability at intraseasonal timescales is much larger than at interannual timescales (Qiu and Chen 2004;Qiu et al. 2009;Cravatte et al. 2015). Sérazin et al. (2015) also show that interannual intrinsic variability is important in this region, suggesting a nonlinear upscaling from mesoscale fluctuations towards lowfrequencies: nonlinear dynamics -including baroclinic instability, scale interactions and rectification (e.g., Sérazin et al, 2018, Zanna et al., 2018 may be strong enough to eventually generate intrinsic interannual variability 85 (Belmadani et al. 2017;Davis et al. 2014). Therefore it is plausible that part of the interannual variability arises from the intrinsic ocean variability alone.
The aim of this paper is to test this latter hypothesis, and to evaluate the respective parts of the external atmospheric variability and of the chaotic oceanic processes in driving the interannual variability of the transports in the Southwest Pacific. For that purpose, a probabilistic modeling strategy is used: an ensemble of 50 global 1/4° 90 ocean/sea-ice simulations, driven by the same realistic atmospheric forcing over 1960-2015 and differing only by slightly perturbed initial conditions (Penduff et al., 2014), is analyzed over 1980-2015 to interpret the currents interannual variability. This so-called "OCCIPUT" oceanic ensemble simulation along with seasonally-forced simulations have been already used to show that internal oceanic dynamics can spontaneously generate a substantial variability over a wide temporal spectrum, including interannual, decadal and long-term trends 95 (Gregorio et al. 2015;Penduff et al. 2018;Serazin et al. 2017;Llovel et al. 2018;Penduff et al. 2019). The oceanic variability that is "deterministically" driven by the prescribed atmospheric variability (e.g. containing the ENSO or SAM atmospheric signal) is estimated from the ensemble mean evolution. Part of the prescribed atmospheric forcing may also originate from chaotic behavior in the coupled system. However as we focus on the ocean, we will use the term "deterministic" to characterize the variability in the ocean common to all 100 simulations. The intrinsic oceanic variability is then quantified from the random dispersion of the 50 members around the ensemble mean. To characterize this variability arising from ocean internal processes, we will use the term "chaotic" in the rest of the paper.
Section 2 describes the ensemble of simulations used, and how deterministic and chaotic variability are quantified. It also describes the data used for validation. Section 3 quantifies the deterministic and chaotic The OCCIPUT ensemble is made of 50 global ocean-sea ice hindcasts, run for 56 years . The model used is NEMO [Madec, 2008], implemented with a ¼° spatial resolution and 75 vertical levels (47 levels 120 in the first 1000m). Model parameters (numerical scheme, subgrid-scale parameterizations) are described in [Bessières et al. 2017]. At the lateral boundaries, a free slip boundary condition is applied.
Each simulation is forced by the realistic Drakkar Forcing Set (DFS) version 5.2, based on the ERA-40 and ERA-Interim reanalyses (Dussin et al., 2016). After a 21-year common spinup , the 50 members of the ensemble are generated in 1960 by activating a small stochastic perturbation in the equation of 125 state within each member [Brankart et al. 2015;Bessières et al. 2017]. This perturbation is only applied for one year: it is switched off at the end of 1960, when the 50 members are restarted from slightly perturbed initial conditions and driven by the same atmospheric forcing. As heat fluxes are computed using bulk formulae, they are slightly different in the 50 members because of sea surface temperature (SST) differences. Wind stresses are computed using absolute winds, without ocean currents feedbacks and are thus identical in each of the members.

130
We focus our analyses on the 1980-2015 period.
The 4D fields of temperature and velocity used in this study are available as monthly means within each member. In addition, temperature and velocity 5-days averages are also available at specific depths for each member. In this study, we focus mainly on two quantities: the 0-1000m integrated transports, and the EKE at the surface and at 100m. The data processing is done as follows: 135 -0-1000m integrated zonal and meridional transports are computed from monthly velocities for each member.
-Long-term trends (very-low frequency signals) are then computed and removed from each time series within each member. This is done using a nonlinear, second-order local regression method (LOESS; Cleveland and Devlin 1988), which high-pass filters time series with a 9-year cutoff period. This method preserves the length of original time series without adverse edge effects (see Serazin et al., 2015) 140 -Detrended transport time series are then low-pass filtered with a 25-months Hanning filter, in order to isolate the interannual variability (such filter eliminates all variability at periods lower than 1 year). Our processing thus confines our analyses and results on time scales between 1 year and 9 years.
We also use the 5-day zonal and meridional velocities at the surface and at 100m depth, filtered at periods lower that 180 days with a Hanning filter, to compute the EKE. EKE monthly anomalies from the mean monthly 145 climatology are then computed, and the long-term trend is removed with a linear trend. Detrended time series are then low-pass filtered with a 25-months Hanning filter, to isolate the interannual EKE modulation.

150
For a given monthly quantity f in member i, we define for each time step t: where * ( ) = * + ∑ " ( ) "-+ "-* is the ensemble mean of the N=50 members (i.e., the atmospherically forced component). " # ( ) is the deviation from the ensemble mean (i.e. the oceanic chaotic component).
For the low-pass filtered signals F in member i, we similarly define for each time step t: The intensity of the monthly forced variability is estimated with the variance of the ensemble mean time series, the bar denoting the time average:  The intensity of the interannual forced variability is estimated with the variance of the ensemble mean time series F * (t) low-ass filtered as explained above:

666666666666666666666666
( 5) where N=50 is the number of ensemble members.
Finally, we define the deterministic variance ratio >1 = 100 * The interannual variances and the RLF ratio depend on the spatial scales considered (Serazin et al., 2015). At a model grid point, they are greater than for a field first averaged over a spatial area (see section 3). Both point-wise and regional diagnoses will be presented; the first diagnosis for quantifying the chaotic variance and deterministic 175 variance ratio in the framework of in situ data comparison, and the second one for more climate-relevant quantities.

Estimation of key current transports
In each member, the 0-1000m transports of key currents are computed as follows (see Figure 1, and Table 1 for a 180 summary) • The transport entering into the Solomon Sea is computed by integrating the 0-1000m meridional transport from coast to coast between 150°E and 162°E, along 10.5°S (labelled 1 on Figure 1a).
• The transport entering into the Coral Sea is computed by integrating the 0-1000m zonal transport along 163°E, from northern tip of New Caledonia to the Solomon Islands (labelled 2 on Figure 1a).

185
• The transport of the NCJ is computed by integrating the 0-1000m westward only zonal transport along 163°E, from 19°S (the northern tip of New Caledonia's reef) to 16.5°S (labelled 3 on Figure 1a); the transport of the SCJ is computed by integrating the 0-1000m westward only zonal transport along 167°E, from 22.5°S (at the southern tip of New Caledonia's reef) to 30°S (labelled 3 on Figure 1a).
• The Tasman Front transport is computed by integrating the 0-1000m eastward only zonal transport 190 along 165°E from 40°S to 30°S (labelled 3 on Figure 1a).
• For the EAC at various latitudes (labelled 5 to 7 on Figure 1a), in addition to the southward flow only, the net transport of the southward flow and the northward adjacent recirculation is computed, following recirculation edge located west of 157.5°E. The total EAC transport is also computed and shown.
Intrinsic variability for both these estimates is also given.

200
To validate the model outputs, we use a mean absolute geostrophic velocity product, based on geostrophic shear estimated from the CARS climatological hydrographic atlas (Ridgway et al. 2002;Condie and Dunn 2006) [http://www.cmar.csiro.au/cars], and a 1000-m absolute velocity field based on Argo float drift (Kessler and Cravatte, 2013a). These 3D zonal and meridional currents are referred to as "Argo-Merged".
We also use the OSCAR (Ocean Surface Current Analysis Real-time) near-surface ocean current estimates 205 (Bonjean and Lagerloef, 2002), from 1993 to 2015. Data are on a 1/3° grid with a 5-day resolution. The horizontal currents are directly estimated from sea surface height, surface vector wind and sea surface temperature. The OSCAR surface current is representative of the top 30m of the upper ocean.
Finally, we use the Nino3.4 index, defined as the average of the SST anomalies in the 5°S-5°N, 170°W-120°W box, from the HadISST dataset (Rayner et al., 2003), and the PDO (Pacific Decadal Oscillation) index. The PDO 210 is one representation of the decadal variability in the Pacific; it exhibits SST spatial patterns similar to those of ENSO in the tropics, but with a broader meridional extent (Mantua and Hare, 2002), and a time series that is primarily decadal. inducing barotropic instability. This weaker EKE may thus be linked to the absence of a CSCC in the model. In other regions, it is worth noting that the model being only eddy-permitting, a significant part of the observed EKE is missing. In fact, mesoscale activity is known to significantly increase in the NEMO model's higher resolution configurations (e.g. Serazin et al., 2015). Finally, the different members exhibit a relatively limited dispersion in 250 mean EKE (dots on Figure 1c, showing regions where the dispersion reaches 10 to 20% of the mean EKE), and have an amplitude similar to that of the ensemble mean EKE.

3.2
Variability of the 0-1000m transport 255 Figure 2a shows the 1980-2015 interannual variance of the 0-1000m zonal transport ensemble mean (i.e, the deterministic, atmospherically-forced, variance of these transports at interannual timescales). Figure 2b shows the intrinsic interannual variance of the 0-1000m zonal transports, and Figure 2c shows the percentage of their deterministic variance. In the ensemble mean, the interannual variability is stronger in the equatorial band, and in the low-latitude western boundary current system, both in the Gulf of Papua and inside the Solomon Sea. In these 260 regions, the variance is clearly deterministic, and atmospheric forcing explains more than 90% of the transports' interannual variance. The forced interannual variability is also strong in the Tasman Front, especially north of New Zealand where it reaches 400 cm 2 .s -2 .
The intrinsic interannual variability of zonal transports is particularly strong in three regions: in the EAC and East Auckland Current (EAUC) system regions, south of 20°S, and east and south of New Caledonia, all along the 265 STCC path. Figure 2c reveals that in these regions, the intrinsic variability strongly dominates the deterministic variability at interannual timescales. This is particularly striking in the EAC extension region, where less than 5% of the interannual variance is deterministic. Clearly, the regions exhibiting strong intrinsic variance coincide with the regions of high EKE, suggesting that mesoscale intrinsic high-frequency variability spontaneously cascades toward lower frequencies ).

270
The same computation is performed after a spatial smoothing of the transports over 10° in longitude and 4° in latitude, to isolate the imprint of interannual intrinsic variability at larger scales ( Figure 2d). This anisotropic filtering takes into account the larger spatial scales in longitude of the zonal currents compared to the latitudinal scales. Although the relative contribution of intrinsic variability tends to be smaller for large-scale signals, the region where it dominates over the deterministic variance remains unchanged. The interannual intrinsic variability 275 https://doi.org/10.5194/os-2020-102 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License.
thus does not only modulate small-scale transports at model grid points, but also has a strong imprint on largescale, climate-relevant transports.

280
Interannual transport anomalies of key currents in the region are also computed in each member (see section 2.1 and Table 1 for details, and Figure  and 7% respectively of these transports interannual variance is chaotic (13%, 9% and 10% for monthly means).
On the other hand, the interannual variability of the transport of the SCJ, the Tasman Front and the EAC at various latitudes is largely chaotic, and strongly impacted by ocean internal processes. At 32°S, near the bifurcation latitude, only 14% of the interannual transport variance is deterministic. In other words, inferring the interannual 290 transport changes of these currents from a single simulation is very uncertain, as it will greatly depend on the initial conditions, and predicting them is very likely to be challenging. Surprisingly, the Tasman Front transport appears as correlated with ENSO, with a 17 months lag. Also surprisingly, at the southern end of the region at 40°S, the southern EAC extension transport is mostly deterministic, although it is known that a substantial fraction of this transport is driven by eddies (van Sebille et al. 2012). Such a feature may arise from a substantial sensitivity of the 295 mesoscale field itself to the interannual atmospheric forcing, which may then rectify the 0-1000m transports in a rather deterministic way. Testing such an hypothesis is left for the future.

Variability of EKE
The same diagnoses are performed for the interannual variability of the EKE at 100m, and are shown in Figure   300 4a,b,c. It is interesting to see that although a part of the EKE modulation at interannual timescales is atmospherically-forced (Figure 4a), a significant part of this interannual modulation is driven by ocean-only internal dynamics (Figure 4b,c). This is especially striking between New Caledonia and New Zealand, in the STCC/Tasman Front regions. These findings are not in agreement with Rieck et al. (2018) who concluded that the EKE low-frequency variability was mostly deterministic in the STCC region. The possible reasons for these 305 disagreements are discussed in Section 6. In the tropical band, the EKE interannual modulation is mostly atmospherically driven.
In addition to interannual variability, the lower-frequency evolution in EKE is also documented. Trends in both

Variability of the EAC separation latitude
Finally, diagnoses are also done for more integrated quantities. The latitude of the EAC bifurcation is computed in each member, following Cetina-Heredia et al. (2014) and Oke et al. (2019a). The SSH isoline contour at 28°S (where the EAC stream is well defined) associated with the maximum southward geostrophic current is identified, and followed southward until it veers eastward. The corresponding latitude is identified as the separation latitude  It confirms that the definition of a bifurcation latitude, and of a well-defined continuous eastward flow emanating from the coast is not an adequate description of the eddying circulation (Oke et al., 2019b). The EKE is also modulated at decadal timescales, in connection with the PDO in the STCC region, as suggested 375 previously by Rieck et al. (2018) and Travis and Qiu (2017) (Figure 7b). During positive phases of the PDO, the EKE decreases in the whole STCC band, in addition to the equatorial band. The processes explaining this decadal deterministic modulation have been investigated in Rieck et al. (2018) and Travis and Qiu (2017). EKE is modulated through a mix of processes including changes in large-scale zonal currents strength, associated vertical shear and stratification. These changes, forced by decadal wind changes in the subtropical gyre, work 380 simultaneously to modify the baroclinic instability intensity, and the EKE levels. As shown in Figure 4c, this deterministic forcing however only represents 10 to 30% of the EKE variability at interannual timescales. It is worth noting that these correlations do not clearly emerge in individual members (see dots in Figure 7).

Drivers of deterministic variability
We further examine the correlation between the vertical shear (defined as the difference between the zonal current at 600m and at 200m) (see Travis and Qiu, 2017), and the EKE levels in the STCC region (160°E-180°, 22°S-385 28°S) (Figure 8). In the ensemble mean, a significant correlation (0.75) is found between the shear and the EKE ( Figure 8c, red dashed line), as also found in altimetric data and a reanalysis in Travis and Qiu (2017) or an OGCM simulation in Rieck et al. (2018). In each individual member however, there is not always a significant correlation between vertical shear and EKE amplitude. The vertical shear appears to be driven at least partly by the atmospheric forcing, as the members do not exhibit a substantial dispersion (Figure 8a), and 80% of the variance 390 is deterministic. On the contrary, the EKE interannual modulation is dominantly intrinsic (Figure 8b), as 40% of the variance only is deterministic. These results highlight the complexity of EKE generation, and the interplay of different forcing processes, whose relative importance varies. The vertical shear, the stratification, and the propagations of eddies from remote regions may all influence the EKE levels. Travis and Qiu (2017) showed that these processes vary in time, and can alternately work to attenuate or reinforce each other. Here, we suggest in 395 https://doi.org/10.5194/os-2020-102 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License. addition that they vary differently in the members, and can alternately attenuate or reinforce each other in the different members.

Spatio-temporal structure of the chaotic oceanic variability
What are the processes at work producing the interannual intrinsic variability? How do mesoscale eddies imprint In all members, the characteristics of the eddies are similar (periods, amplitude), as expected since they share the same model dynamics. However, their phasing differs, and their interannual modulation phase also differs from one member to the other. Further investigating the non-linear interactions leading to the transport interannual variability, probably arising from eddy low-frequency rectification, would be interesting but is beyond 420 the scope of this study.
In the STCC region, the spatial structure of the chaotic variability is quite different (Figure 10 a, b). In terms of zonal transport, it consists in zonally elongated structures, alternating meridionally, slanted in a southwestnortheast direction, and slowly propagating southward, as can be seen when comparing the two first EOFs and as weak latent jets that shows up only on filtered signals such as the EOF decomposition performed here.
Interestingly, monthly and interannual EOFs and PCs are similar (Figure 10c), revealing that large-scale intrinsic variability in this area is dominated by a low frequency (3-4 year timescale) signal. These anomalies are not in phase in the different members, confirming that they arise from ocean-driven intrinsic processes, potentially through nonlinear scale interactions (i.e., an inverse cascade) as discussed above.

440
In this study, the origin and features of interannual variability of the transports and EKE in the Southwest Pacific have been investigated using a probabilistic modeling strategy. A large ensemble of 50 simulations performed with the same OGCM at eddy permitting resolution, perturbed initially but driven by the same realistic atmospheric forcing over 1980-2015, has been analyzed to disentangle the atmospherically-forced, deterministic ocean variability and the internal, chaotic ocean variability. This approach had previously been used to show that 445 chaotic oceanic processes can have more impact than the deterministic atmospheric variability on the lowfrequency variability and trends of regional transport, sea level or oceanic heat content, particularly at mid-latitudes The spatial-temporal features of the oceanic chaotic variability are complex, and varies from region to region. Interestingly, as the simulations share the same dynamics, they produce chaotic structures with similar spatial characteristics but with different phasing. In the STCC, the low-frequency imprint of these structures is zonally elongated, and its spatial organization is likely the result of an incomplete anisotropic inverse cascade due 460 to nonlinear scale interactions (e.g., Sérazin et al., 2018). In other regions, they are eddy-shaped. The physical processes governing this chaotic variability remain to be fully understood.
Our results on EKE contrast with the findings of Rieck et al. (2018), who concluded that the EKE decadal variability in the STCC region was, uniquely among the subtropical gyres, mostly deterministic, driven by windforced decadal changes. In agreement with Travis and Qiu (2017) analyses of observations, they explained these 465 decadal EKE modulation by decadal changes in vertical shear, driven by wind stress curl changes in the South Pacific, while still recognizing that other processes such as changes in stratification or density anomalies remotely forced also contribute to modulate the EKE low-frequency changes. The disagreement between our and these authors conclusions is likely due to our different modeling strategy, which gives access to both intrinsic and EKE generation, and the interplay of different forcing processes; they also emphasize the importance of taking simultaneously into account the chaotic ocean behavior and the imprint of the atmospheric variability in an ensemble approach. It is also worth noting that our results do not question the existence of deterministic, forced interannual variability of EKE in the area; our results mostly emphasize that the chaotic ocean interannual variability may overwhelm this deterministic variability over large regions.

480
The fraction of interannual variability explained by chaotic oceanic processes depends on the spatial scale considered. At model grid points, this fraction is greater than for spatially integrated variables, or even for largescale transports across sections. Yet, our results show that for some key currents, section-integrated transport variability is dominantly chaotic. Moreover, our study reveals that chaotic variability may result in large-scale, spatially coherent signals at interannual timescales. The EAC separation latitude variability (defined where eddies 485 detach off the EAC) at interannual timescales is also largely chaotic.
These results have important implications for the understanding and forecasting of the transports in the region, and for observational strategies. Firstly, they show that the impact of climate modes of variability such as ENSO or PDO on the interannual variability of transports and EKE may be hidden by oceanic internal variability, and that an ensemble of simulations is needed to bring out this signature. This clearly hampers our capacity to 490 attribute the forced origin of the interannual variability of transports from a single simulation.
Importantly, these results highlight that interpretation of in situ or satellite observations in areas of high intrinsic variance should be done with caution. While the low-frequency signals amplitude does not strongly vary between members, their phase may be highly random because of chaotic oceanic processes, and may be largely uncorrelated with the atmospheric variability. Figure  It should be kept in mind that we used an eddy-permitting (1/4°) resolution that sits at the limit of presently affordable long-term global ocean large ensembles, but in which only part of the mesoscale activity is simulated.

510
It is probable that at higher resolution, with the increased mesoscale activity, inverse cascades of energy toward large scales and low frequencies are stronger: this has been shown in terms of sea level by Serazin et al. (2015), who suggested a direct feeding of low-frequency chaotic variability by mesoscale activity both within and away from eddy-active regions. Gregorio et al (2015), however, have suggested that the forced variability of the AMOC is also likely to increase when switching from 1/4° to 1/12° resolution; whether this also holds for the variables 515 https://doi.org/10.5194/os-2020-102 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License.
investigated here needs to be verified, but these authors did not report a substantial sensitivity of forced-to-total variability ratios (noted RLF here) to this resolution increase. More generally, it would be interesting to assess whether the OCCIPUT ensemble gives a robust estimate of the ocean chaotic variability; comparisons with other ensemble simulations (such as that performed at 1/10° by Nonaka et al. (2020) would be of use for such an assessment.

520
We also used forced and not coupled oceanic simulations, which is an indispensable step to isolate the chaotic ocean-only contributions. The importance of air-sea interactions at mesoscale are however known to be significant.
The wind stress variability is impacted by the ocean current feedback, which could damp the EKE eddy-scale variability in return (e.g. Renault et al. 2017Renault et al. , 2019. Intrinsic variability is also likely to impact the air-sea heat flux variability through sea surface temperature fluctuations and affect the atmospheric variability in return. These 525 possible effects should be investigated in the future as the observed variability reflects the fully coupled climate system. https://doi.org/10.5194/os-2020-102 Preprint. Discussion started: 16 October 2020 c Author(s) 2020. CC BY 4.0 License.
Data availability: The model dataset used for this study is available on request (thierry.Penduff@cnrs.fr).
Author contributions: Thierry Penduff lead the OCCIPUT project and run the experiments. Sophie Cravatte performed the analyses, prepared the manuscript, with contributions from all co-authors. All authors have contributed to the interpretations of the results and to a significant part to the presented scientific work.
Competing interests: The authors declare that they have no conflict of interest.