## Abstract

An intrinsic mode of self-sustained, interannual variability is identified in a coarse-resolution ocean model forced by an annually repeating atmospheric state. The variability has maximum loading in the Indian Ocean, with a significant projection into the South Atlantic Ocean. It is argued that this intrinsic mode is caused by baroclinic instability of the model’s Leeuwin Current, which radiates out to the tropical Indian and South Atlantic Oceans as long Rossby waves at a period of 4 yr. This previously undescribed mode has a remarkably narrowband time series. However, the variability is not synchronized with the annual cycle; the phase of the oscillation varies chaotically on decadal time scales. The presence of this internal mode reduces the predictability of the ocean circulation by obscuring the response to forcing or initial condition perturbations. The signature of this mode can be seen in higher-resolution global ocean models driven by high-frequency atmospheric forcing, but altimeter and assimilation analyses do not show obvious signatures of such a mode, perhaps because of insufficient duration.

## 1. Introduction

The ocean exhibits variability on a range of space and time scales, from centimeters to planetary scales, and from seconds to millennia. Variability on micro- and mesoscales (hundreds of kilometers and less) is essentially turbulent and unpredictable in detail (Wunsch 2001). Variability at lower frequencies and larger spatial scales is often more organized and can take the form of recognizable modes or oscillations. Familiar examples of such modes include El Niño–Southern Oscillation and the Atlantic multidecadal oscillation as well as a host of gyre, basin, and thermohaline oscillations (Dijkstra and Ghil 2005). Such modes can be sustained by atmospheric noise or air–sea coupling (e.g., Sévellec et al. 2009; Frankcombe et al. 2009) or can take the form of self-sustained intrinsic variability, independent of variability in the atmospheric forcing (e.g., Dijkstra and Ghil 2005; Simonnet et al. 2005; Penduff et al. 2011).

The variability of the ocean circulation due to large-scale modes has implications for the predictability of the ocean and—through its coupling to the atmosphere—the global climate. Large-scale modes have potentially reproducible patterns, possibly enhancing the predictability of the forced response. On the other hand, prediction of the phase at which the pattern will be excited by atmospheric forcing is challenging; this phase uncertainty degrades the predictability of the forced response, which becomes masked by the intrinsic variability.

Intrinsic, low-frequency variability in sea level anomalies (SLA) has been documented in a global, eddy-resolving model (Sérazin et al. 2015) forced by a repeated, climatological atmospheric seasonal cycle. Interestingly, even the large-scale, low-frequency intrinsic variability in this configuration is largely concentrated in regions of high eddy activity, where the small-scale, low-frequency fluctuations are also dominant: these are the western boundary current extensions and the Antarctic Circumpolar Current. In addition to these notoriously eddy-active regions, the southern Indian Ocean exhibits enhanced low-frequency variability, which is eastern, rather than western, intensified.

In this work, we show that Indian Ocean intrinsic variability with the same characteristics as that found by Sérazin et al. (2015) is obtained in similarly forced, low-resolution (1°), global computations. At low resolution, the main source of intrinsic interannual variability is a near-global mode originating from a radiating instability of the model Leeuwin Current (LC), the eastern boundary current of Australia, with an expression in the Atlantic but not the Pacific. Like the baroclinic instability of a western boundary current extension, the primary instability of this eastern boundary current is a vertical shear instability, analogous to that described in Hristova et al. (2008), with the energy source being the available potential energy of the time-mean LC. Because the current is on the eastern boundary, it can radiate large-scale Rossby waves at interannual periods, as long as perturbations originating from the eastern boundary satisfy the free Rossby wave dispersion relation in the interior. While this type of radiating instability can, in principle, originate from any eastern boundary current, we argue that the LC is preferred because of its strength, persistence during all seasons, and coherence along the coast of Western Australia.

The formulation of the model is described in section 2. We argue in section 3 that the dominant mode of variability is organized into at least two beams of baroclinic Rossby waves with a 4-yr period, propagating westward across the south Indian Ocean and into the Atlantic Ocean. The mechanism responsible for these Rossby waves—radiating baroclinic instability of the LC—is identified in section 4 and illustrated using a 2.5-layer quasigeostrophic model. The implications for oceanic predictability as well as a survey of evidence for (and against) the existence of the 4-yr mode is given in section 5. A short summary is provided in section 6.

## 2. The model

The ocean model—referred to here as latitude–longitude–polar cap (LLC)–CORE normal-year forcing (CNYF)—is based on the ECCOv4 (Forget et al. 2015) configuration of the MIT general circulation model (MITgcm; Marshall et al. 1997a,b). Specifically, this is a global configuration with an LLC topology, a 1° nominal horizontal resolution that telescopes to ⅓° near the equator and 50 vertical levels. The physics include dynamic sea ice and the KPP mixed layer parameterization (Large et al. 1994). The Gent–McWilliams (GM; Gent and McWilliams 1990) parameterization of mesoscale eddies is implemented using the boundary value problem approach of Ferrari et al. (2010). We use spatially constant GM and isopycnal (Redi 1982) diffusivities of 500 m^{2} s^{−1}.^{1} In the interior, the diapycnal diffusivity is quite small (*κ*_{υ} = 10^{−5} m^{2} s^{−1}), and the advection scheme is the nearly adiabatic scheme of Prather (1986). The surface forcing is obtained from the annually repeating CNYF atmospheric state of Large and Yeager (2009), version 2, using the bulk formulae of Large and Yeager (2004). Additional forcing is provided by restoring sea surface salinity (SSS) to an annually repeating monthly SSS climatology described in Forget et al. (2015), over a time scale of 300 days. This SSS relaxation rate is intermediate between the moderate and weak salinity restoring cases described in the CORE protocol (Griffies et al. 2009). The total forcing thus contains a broad spectrum of variability whose longest period is 1 yr.

The LLC–CNYF model was spun up to a statistical steady state over approximately 4800 yr. At the end of the spinup period, the maximum trends were less than 2 × 10^{−3} °C century^{−1} for temperature (at 1500 m) and less than 4 × 10^{−4} psu century^{−1} for salinity (at 350 m). Further, the secular variation in transport of transport of all major current systems (the ACC, AMOC, Florida Current, etc.) was less than 0.06% century^{−1}. The analysis in this paper makes use of the last 300 yr of the simulation.

## 3. Intrinsic variability

### a. Characteristics

The interannual component *x*′ of a given field *x* is isolated by removing its monthly climatology. Since the model’s forcing repeats annually, any variability with periods longer than a year is intrinsic to the ocean dynamics. The Indian Ocean stands out as a region of enhanced interannual variability of sea surface height (SSH), with two elongated lobes of high variability (Fig. 1a). The northern maximum in variability extends from the southwestern tip of Australia in a curving arc to approximately 20°S and has its greatest expression in the eastern and central Indian Ocean. In contrast, the southern variability maximum found near 25°S is oriented zonally and is largest in the western Indian Ocean. This variability is not strictly a surface phenomenon; other dynamical quantities such as potential temperature, salinity, and velocity have similar patterns of variability in the upper 1000 m of the water column. The fingerprint of the intrinsic, interannual variability is particularly clear in potential density referenced to 2000 dbar: (Fig. 1b). The North and South Atlantic Oceans also have localized regions of enhanced interannual variability, though of a smaller magnitude than the Indian Ocean. Without interannual variability in the atmospheric forcing, the Pacific basin is quiescent and displays little intrinsic variability.

The two lobes of enhanced variability in the Indian Ocean are associated with two beams of westward-propagating disturbances with maximum amplitudes at distinct latitudes and depths (see animation 1 in the supplemental material). Though they both originate off Western Australia and are temporally synchronized, the two beams occupy distinct depth and density ranges (Fig. 2). Both beams have multiple maxima in the vertical separated by regions of low variability, indicative of baroclinic modal structure. Both beams have significant contributions from the second baroclinic mode, and the southern beam appears to have an additional contribution from the third baroclinic mode. The northern beam has its maximum amplitude at 150-m depth, and the southern beam is most pronounced at 300 m; the beams are less distinctly separated below 400-m depth. (The maximum of the southern beam is found near 200 m in Fig. 2 but becomes deeper toward the west.) As illustrated in Fig. 2, the regions of enhanced variability associated with both beams coincide with local minima in the large-scale Ertel potential vorticity (EPV), defined as

where *N*^{2} is the Brunt–Väisälä frequency. At fixed latitude, Π is proportional to *N*^{2}, so Fig. 2 indicates that the variability in is largest where the background stratification is strongest (*f* is negative in this region). This is consistent with the vertical structure of baroclinic density modes, which attain their maximum amplitude where *N*^{2} is largest (e.g., Ferrari et al. 2010). In the upper water column, the two beams are separated by a pool of near-zero EPV water. On the south side of the pool, the meridional EPV gradient is reversed relative to the planetary potential vorticity gradient. This reversed EPV gradient impedes westward propagation and could be the reason the variability separates into two beams.

Frequency–zonal wavenumber spectra for the boxes and depths shown in Figs. 3a and 3b are calculated by Fourier transforming in longitude and time after weighting by a Tukey window with a flat, central portion spanning 50% of the data. The resulting power spectra, shown in Figs. 3c and 3d, are averaged meridionally over the boxes. Both beams have strong peaks with a 4-yr period; the northern beam has a secondary peak with a period of 2 yr. The peak at 4 yr is not a pure spectral line; while the peak is centered on 4 yr, it has significant^{2} power in the period band 3.8–4.2 yr with a local minimum at a period of exactly 4 yr (visible as a white stripe through the peaks in Figs. 3c and 3d). Spectral broadening of this type can result from a purely harmonic carrier signal with a slowly modulated amplitude and phase; the local minimum at 4 yr implies that the modulating signal is approximately bimodal.

The wavenumber spectra in Figs. 3c and 3d show that the northern and southern beams have wavelengths of approximately 2500 and 2100 km, respectively, with some zonal variation. Comparison with the linear dispersion relation for long Rossby waves

where *L*_{n} is the *n*th baroclinic Rossby radius,^{3} shows that the frequencies and wavelengths are consistent with baroclinic mode-2 Rossby waves, with small contributions from baroclinic mode 3. The waves are shorter and slower at higher latitudes, consistent with the variation of the internal Rossby radii with latitude.

Figure 4 shows the vertically integrated variance of *σ*_{2} obtained via complex demodulation (e.g., Bloomfield 2004, his section 8). The demodulation is performed using a carrier signal with a period of 4 yr and zero phase at 1 January and low-pass filtered using a twelfth-order Butterworth filter with a cutoff period of 8 yr. The first and last 10 yr of the time series are discarded after filtering in order to avoid edge effects. The spatial distribution of potential density variance after demodulation is qualitatively similar to that of the full interannual field (cf. Fig. 1b); the demodulated field explains 78% of the total variance of . We therefore attribute the majority of the interannual variance of to an oscillation with a 4-yr period whose amplitude and phase are modulated on longer time scales.

### b. CEOF analysis

The observed interannual variability is efficiently described by the leading three-dimensional complex empirical orthogonal function (CEOF) **e**_{1}. This CEOF is the leading eigenvector of the complex covariance matrix of , where is complexified using the frequency domain representation of the Hilbert transform (e.g., von Storch and Zwiers 2001, their section 16.3). The discrete values of taken from the model are weighted by the volume of the grid cell in which they reside, so that **e**_{1} optimizes the globally integrated variance of .

The CEOF **e**_{1} captures 35% of the global variance. As shown in Fig. 5, this CEOF is concentrated in the Indian Ocean (where it explains 42% of the local variance), with a significant projection into the South Atlantic Ocean (explaining 31% of the local variance). The amplitude in the Pacific is negligible. The leading CEOF also describes some variability in the North Atlantic Ocean, although the total North Atlantic variability explained by this CEOF is small (<2%) and the physical connection between the Southern and Northern Hemisphere expressions of this CEOF is not clear. In the following, we concentrate on the Indian and South Atlantic Oceans.

The principal component (PC) of this CEOF, *p*_{1}(*t*), is an almost sinusoidal oscillation with a 4-yr period (Fig. 6a). The power spectrum of *p*_{1} is similar to that shown in Fig. 3; it is strongly peaked near a period of 4 yr, but the peak is spread over several frequency bins with a local minimum at a period of precisely 4 yr. The real and imaginary parts of *p*_{1}—*p*_{r1} and *p*_{i1}, respectively—are 90° out of phase on average, although the precise phase of *p*_{1} varies like a Gaussian centered on 90° with a standard deviation of 15°. Since *p*_{r1} and *p*_{i1} are in approximate quadrature, **e**_{1} describes a propagating disturbance. Complex demodulation (as described in section 3a) shows that the amplitude and phase of *p*_{1} varies slowly on decadal and longer time scales (Figs. 6b,c); this slowly modulated 4-yr signal explains 90% of the variance of *p*_{1}.

The nearly exact 4-yr periodicity suggests that the variability might be synchronized with the annual cycle, but—if synchronization were perfect—the extrema of *p*_{1} would occur at the same time of year. However, this is not the case; while 50% of *p*_{r1}’s maxima occur in spring (March–May), these periods of relatively stable phase last only for 2–3 oscillations (8–12 yr) and are separated by intervals when the phase drifts either forward or backward over the course of 1–2 oscillations (4–8 yr; Fig. 6c). The unfiltered intrinsic variability displays a similar lack of seasonal synchronization as the CEOF; wherever the intrinsic variability is large, its phase varies in a manner consistent with the first PC. The 4-yr period therefore appears to be coincidental and not a result of coupling to the annual cycle.

There are two episodes—centered at *t* ≈ 133 and 234 yr—where the amplitude of *p*_{1} dramatically decreases and the phase shifts rapidly over the course of 10–20 yr (Fig. 6). The phase shifts by ~180° during the first episode and ~90° during the second episode. There is an additional hint of a similar, but less dramatic, episode near *t* ≈ 50 yr. The phase is ~60° for almost half of the time series and −120° for the other half, so the phase distribution is approximately bimodal. This bimodalilty of the PC is shared with the full time series and explains why the power spectrum shown in Fig. 3 has peaks on either side of the 4-yr period but a local minimum at exactly 4 yr.^{4}

During the low-amplitude episodes, the second PC *p*_{2}, which explains only 6.4% of the global variance, increases in magnitude and develops a 4-yr oscillation not seen during other times (not shown). The corresponding CEOF **e**_{2} is concentrated in the Indian Ocean with a spatial pattern that is visually similar (though orthogonal) to **e**_{1}. In fact, all of the first 10 CEOFs have significant amplitude in the Indian Ocean, implying that the intrinsic mode may be too complex to be completely described by a single CEOF.

Variations in the phase and amplitude of *p*_{1} appear to display some periodicity on the time scales of approximately 100 yr—roughly the period between low-amplitude episodes—but the periodicity is not exact. The presence of a centennial-scale oscillation cannot be confirmed without considering a longer optimization period for the CEOFs, since Hilbert transforms become contaminated by edge effects when applied to signals with periods comparable to the length of the time series.

In the Indian Ocean, the physical structure of **e**_{1} corresponds to the Rossby wave beams described in section 3a. The CEOF calculation’s inherent signal enhancement clarifies the beam structure compared to the raw data. In particular, it is possible to approximately trace the beam paths by determining, at a given depth, the latitude of maximum amplitude of |**e**_{1}| for each longitude. In general, a ray path is perpendicular to lines of constant phase; this property was verified by visual inspection and minor adjustments made by hand in cases where the originally determined path was not orthogonal to phase lines. Two primary beam paths were identified in this manner and are discussed below.

#### 1) Interbasin beam

Figure 7 shows the spatial structure of **e**_{1} at 300-m depth (see also supplemental animation 2). The real part **e**_{r1} is a train of disturbances with alternating sign with the largest magnitude near 25°S in the Indian Ocean. The extrema of the imaginary part **e**_{i1} are located between the extrema of **e**_{r1}, with the maxima of **e**_{i1} to the east of the maxima of **e**_{r1}. Since *p*_{i1} lags *p*_{r1}, the arrangement of **e**_{r1} and **e**_{i1} indicates westward propagation consistent with the Rossby wave interpretation advanced in section 3a. The train of disturbances originates off the southwest coast of Australia and propagates westward, arcing toward Madagascar. Lines of constant phase are oriented northeast–southwest in the eastern part of the basin and rotate counterclockwise as the wave propagates west, becoming oriented northwest–southeast near Madagascar. These features are consistent with a long Rossby wave undergoing beta refraction.

Upon reaching Madagascar, the disturbances become entrained in the Agulhas Current and are advected around the southern tip of Africa into the South Atlantic Ocean, where they organize into another Rossby wave–like train of disturbances with a smaller amplitude and longer wavelength than the original disturbances in the Indian Ocean. The disturbances propagate northwestward across the South Atlantic Ocean, eventually reaching the coast of South America, where they convert into coastally trapped waves and propagate toward the equator. The disturbances may influence variability in the North Atlantic, but the causal connection between variability in the North and South Atlantic in this model is unclear. The vertical structure of **e**_{1} is a mixture of the second and third baroclinic modes, with the third baroclinic mode more prominent in the Indian Ocean than the South Atlantic. The modal structure is visually apparent in Fig. 8 and is confirmed by explicit decomposition of the disturbances into vertical modes (not shown). The disturbances are much shallower in the South Atlantic Ocean than in the Indian Ocean, with the extrema below 600 m entirely absent. This loss of third-mode structure may be simply because high baroclinic modes are damped more rapidly than low modes by viscosity and diffusion.

Wavelengths and phase speeds can be estimated from the distances between extrema and the known 4-yr periodicity. In general, the waves are slower at higher latitudes (approximately 0.6 cm s^{−1} at 20°S and 0.2 cm s^{−1} at 30°S), consistent with the latitudinal dependence of the long Rossby wave phase speed. The wavelength of the disturbances varies proportionally with the phase speed to maintain a constant period, ranging from 1200 km near Australia to 2500 km east of Madagascar to 3400 km in the lowest latitudes of the South Atlantic. This is consistent with the frequency–wavenumber spectra discussed in section 3a. The propagation time scale of this mode of intrinsic variability can be estimated by counting the nodes of the CEOF and multiplying by the period. This gives about 32 yr from the southeastern Indian Ocean off Western Australia to the equatorial Atlantic.

#### 2) Northern beam

The second most prominent train of disturbances is found north of those discussed in the previous section and is concentrated near 150-m depth. This northern beam also consists of westward-propagating disturbances with the characteristics of Rossby waves; these waves appear to originate near the southwestern tip of Australia, propagate in a northwestern arc across the Indian Ocean, and eventually encounter the northern tip of Madagascar (Fig. 9). The vertical structure of this Rossby wave train is dominated by the second baroclinic mode (Fig. 10). As with the interbasin beam, the zonal phase speed of this beam increases as the disturbances propagate into low latitudes, while the wavelength increases to keep the period constant. Since this beam propagates farther north in the Indian basin than the interbasin beam, the increase in zonal scales is more dramatic; the wavelength grows from approximately 400 km near the southwestern tip of Australia to more than 3500 km when the disturbances encounter Madagascar.

Once these disturbances reach Madagascar, they refract around its northern tip, radiate toward Africa, and propagate as coastally trapped waves up the coast of Africa to the equator. The waves then quickly cross the basin to Malaysia in the equatorial waveguide (likely as equatorial Kelvin waves), as shown in Fig. 11. These features suggest that the northern beam might be part of a basin mode concentrated in the Indian Ocean.

#### 3) Beam steering

The discussion in the preceding two sections is framed in terms of Rossby waves propagating through a stationary medium. Of course, the southern Indian Ocean is host to a horizontally and vertically sheared wind-driven circulation (e.g., Talley et al. 2011, their section 11). Determining the precise impact of the wind-driven background flow on low-mode baroclinic Rossby waves is not straightforward; geostrophically balanced currents alter the background potential vorticity gradient, so a Rossby wave propagating in a background flow is not simply Doppler shifted by the mean current (this is the non-Doppler effect).

Qualitatively, the Rossby wave beams align with the streamlines of the horizontal velocity field averaged over the deepest maxima in their amplitudes: 275–600 m for the northern beam and 600–950 m for the interbasin beam (Fig. 12). The mean velocities at these depths are westward and comparable to the phase speeds of the Rossby waves themselves. Advection by these mean currents accelerates the Rossby waves, and the wavenumber and amplitude decrease to conserve period and wave action (Andrews and McIntyre 1978). This effect is especially pronounced in the northern beam (cf. Figs. 9 and 10) where the increase in wavelength exceeds that which would be expected from its equatorward propagation alone. In the western Indian Ocean, the phase speeds of the waves are intermediate between the mean currents shown in Fig. 12 and the intrinsic Rossby wave phase speeds due to a combination of the Doppler and non-Doppler effects.

## 4. Mechanism: Baroclinic instability of an eastern boundary current

Insight into the mechanism responsible for this intrinsic mode of variability can be found by considering the eddy kinetic energy (EKE) budget, which is derived in appendix A. The two significant source terms in the EKE budget are the baroclinic and barotropic conversion rates

and

respectively. The overbar indicates the time mean, a prime indicates the deviation from the time mean, **u** is the horizontal velocity vector, is the unit vector in the vertical direction, *ϑ* is the latitude, and *a* is Earth’s radius.

The geographical distribution of the two conversion terms is shown in Fig. 13. Conversion of the mean to eddy kinetic energy is large primarily near the equator, where the refinement in grid resolution allows the formation of narrow currents susceptible to horizontal shear instabilities. The conversion of potential energy into eddy kinetic energy is also enhanced near the equator and large near the major western boundary currents. In addition, a prominent local maximum of is associated with the LC west of Australia, and is generally positive throughout the south Indian Ocean. It thus appears that the intrinsic variability observed in the south Indian Ocean is sustained by conversion of potential energy into eddy kinetic energy.

The mechanism leading to the spontaneous excitation of the intrinsic mode is, we believe, a radiating instability of the model LC off Western Australia. The instability mechanism is baroclinic, even though the radius of deformation is not resolved in this model. Both baroclinic and barotropic radiating instabilities can be generated by eastern boundary currents (Hristova et al. 2008; Wang et al. 2013), and in our case the baroclinic conversion rate (Fig. 13a) in the Indian Ocean is large, while the barotropic conversion rate is negligible (Fig. 13b). Although other eastern boundary currents are present, the model LC is strong, persistent, and unidirectional along the entirety of the west Australian coast. Unlike other boundary currents whose surface expression changes sign seasonally, the surface LC is poleward year-round and has an equally strong equatorward Leeuwin Undercurrent (LUC).

Radiating baroclinic instabilities of eastern boundary currents have been discussed by Hristova et al. (2008); a vertically sheared, purely baroclinic, steady, meridional eastern boundary current is unstable to radiating waves whose wavelengths in the alongshore direction exceed the Rossby deformation radius. These unstable modes radiate into the interior as long Rossby waves with a westward group velocity.

Because of the presence of the LUC, the time-mean eastern boundary current has a baroclinic structure with at least two nodes in the vertical (baroclinic mode 2). A possible idealization of this configuration is to consider the stability of a current with a single baroclinic-mode vertical structure, subject to perturbations with a vertical structure with two baroclinic modes. This system can be modeled using a 2.5-layer quasigeostrophic (QG) model, with a prescribed steady meridional current confined to the eastern boundary of a semi-infinite domain. Because we are interested in disturbances spontaneously generated on the eastern boundary and freely propagating westward, we neglect the influence of interior mean flows that could alter the stability and propagation properties of the unstable modes.

The details of the linear stability analysis are given in appendix B. A systematic study of the 2.5-layer instability of an eastern boundary current is beyond the scope of this work. Here, we simply note that the unstable solutions (cf. Fig. 14) have low frequencies (period about 4 yr), relatively long meridional wavelengths (400 km), and radiate toward the west as long Rossby waves (zonal wavelengths ≈1800 km). The properties of these waves closely match those of the dominant waves observed in the LLC–CNYF model (cf. Fig. 3). For the unstable modes, the frequency of the radiated Rossby waves will have a positive imaginary part and the radiated long Rossby waves will decay in space toward the west, while the neutral modes will maintain their amplitude. Presumably the nonlinearly equilibrated waves also maintain their amplitude as the Rossby waves travel westward. The essential point is that the radiated waves have scales much larger than the deformation radii, even though this scale is generally implicated in baroclinic instability. (The unstable solution shown in Fig. 14 has some structure on the scale of the Rossby radius within the localized eastern boundary current, but the radiated waves have much larger scales.)

In addition to the low-frequency modes directly excited by baroclinic instability, it is also possible that the mechanism of subharmonic excitation of secondary waves—discussed by Wang et al. (2013) in the context of barotropic dynamics—is operative in the baroclinic case. In the barotropic case, a vertically uniform, steady, meridional eastern boundary current with horizontal shear is unstable to modes that undergo a secondary instability to waves with double the period and wavelength that radiate into the interior. In the barotropic context, the period of the radiated Rossby waves is appropriate for free barotropic Rossby waves and thus much shorter than in the baroclinic case.

This linearized model is meridionally periodic with a meridionally uniform background state, resulting in the generation of plane waves for the instability. The primitive equation model run shows at least two separated beams, or paths, that follow different routes across the Indian Ocean. This collimation is likely due to inhomogeneities of the Indian Ocean waveguide—particularly, the pool of near-zero EPV water shown in Fig. 2—which is not represented by the QG model.

## 5. Discussion

### a. Implications for decadal predictability

Because of its large heat capacity and relatively slow evolution, the ocean has a longer memory than the atmosphere. If the ocean is itself predictable, its long memory may be used to enhance predictions of regional climate variation on interannual to decadal time scales.

The effect of intrinsic variability on oceanic predictability depends on the nature of the variability. A regular oscillation would enhance predictability since the future evolution of the oscillation is determined once its phase is. Irregular variability, on the other hand, reduces predictability by obscuring the ocean’s response to changes in forcing or initial conditions. While the intrinsic variability described in this paper has aspects that appear regular (it is highly peaked at a 4-yr period), its amplitude and phase vary chaotically on decadal time scales. A generic feature of chaotic systems is that nearby trajectories diverge exponentially for small separations. Chaotic intrinsic variability thus appears as an instability of the time-varying background flow. More precisely, if the original model trajectory is denoted by **x**_{0}(*t*) and **x**_{1}(*t*) is a nearby trajectory that differs from **x**_{0} by a small perturbation in initial conditions or forcing, then the difference

will initially grow exponentially if the system is chaotic. The growth of *δ***x** will eventually saturate when the separation between the two trajectories is comparable to the size of the system’s attractor.

To test whether trajectories diverge exponentially, the model was restarted using pickup files that were truncated from double- to single-precision floating point numbers. This truncation is equivalent to a random perturbation of relative magnitude 10^{−7} applied uniformly to all variables. The blue line in Fig. 15 shows the magnitude of *δ***x** as measured by the squared difference in potential density *σ*_{2}, averaged over the global ocean. After an initial transient, the *δ***x** grows exponentially with an *e*-folding time of 3 yr over a decadal time scale, eventually saturating after approximately 20 yr. The difference *δ***x** has a strong projection on CEOF 1 (green line in Fig. 15).^{5} In fact, the projection onto CEOF 1 grows at a faster rate than *δ***x** itself. The projection saturates at about one-third the magnitude of *δ***x**, which is consistent with the variance of the background state explained by CEOF 1. These results indicate that the ocean model is, indeed, chaotic and that much of the chaotic evolution is driven by the intrinsic variability centered in the Indian Ocean.

The extent to which the intrinsic mode influences regional climate variations depends on whether it strongly influences the exchange of properties and energy—particularly heat—between the ocean and atmosphere. Snapshots of the surface heat flux anomaly *δQ* = *Q*(**x**_{1}) − *Q*(**x**_{0}) for the initial condition perturbation experiment are shown in the left panels of Fig. 16. The values of *δQ* remain small in the Pacific, so only the Atlantic and Indian Oceans are shown. Each panel shows August, which is near the end of austral winter when surface heat flux anomalies are largest in the Southern Hemisphere. The first panel is 8 months after the initial perturbation; *δQ* is very small and concentrated near the equator, the Atlantic western boundary currents, and the southwestern corner of Australia. By year four, *δQ* is larger by a factor of 50 and is beginning to form alternating patches of positive and negative anomalies in the southeastern Indian Ocean. These anomalies propagate westward and amplify over the next decade, crossing most of the Indian Ocean by year 10. The lower-left panel of Fig. 16 shows year 30, by which time the instability has saturated with heat flux anomalies exceeding 30 W m^{−2}. The saturated anomaly is statistically indistinguishable from an interannual heat flux anomaly taken from a random August of the unperturbed model. Note that most of the heat flux anomaly is concentrated in the midlatitudes and falls to negligible values (<2 W m^{−2}) by 20°S; thermal anomalies due to the intrinsic mode are thus expected to primarily affect the climate of the midlatitudes (e.g., Australia and Madagascar) rather than the tropics.

The intrinsic variability also dominates the response to forced perturbations, as shown in the right panels of Fig. 16. In this experiment, a heat flux perturbation with a maximum amplitude −30 W m^{−2} is applied over the Labrador Sea over the course of a randomly chosen winter season. The spatial footprint of the perturbation is a Gaussian with a half-width at half maximum (HWHM) of 450 km centered on 60°N, 56°W, and its temporal structure is given by a Gaussian centered on 15 February with a HWHM of 45.5 days. The integrated thermal anomaly associated with this perturbation is 2 × 10^{15} J. The local response in the Labrador Sea is visible in the upper-right panel but weakens and fragments into an incoherent patch of alternating anomalies in the North Atlantic by year four. The Indian Ocean response, starting in the LC off Western Australia, is apparent by year 10—too fast to be affected by advection or propagation of baroclinic waves from the Labrador Sea. The rapid response of the Indian Ocean is likely due to the nonlocal nature of the elliptic equation used to determine SSH, which causes changes in SSH to be felt everywhere instantaneously. The change in SSH far from the site of the perturbation will be small—on the order of numerical roundoff error—if the elliptic solver is accurate. However, since the system is chaotic, these small changes will be amplified exponentially and affect the future evolution of the system. Indeed, after 30 yr, the saturated *δQ* pattern is statistically indistinguishable from that resulting from an initial condition perturbation (cf. the bottom panels of Fig. 16).

These two examples show that the intrinsic variability found in this model is sufficient to overwhelm the response to a forcing perturbation on decadal time scales. Unless the detailed temporal evolution of the intrinsic mode can be understood and predicted, its presence severely limits oceanic predictability. Further, as models used for climate prediction move to higher resolutions and become less dissipative, similar large-scale intrinsic variability could affect the detection and attribution of anthropogenic climate change by masking the forced response of these models (e.g., Richter and Marzeion 2014; Lyu et al. 2015; Sérazin et al. 2016). The fact that vigorous intrinsic variability is observed in the LLC–CNYF model, which has a similar resolution but is less dissipative than a typical climate forecast model, suggests that intrinsic oceanic variability may become significant in climate models long before they attain eddy-resolving resolutions.

### b. Comparison with observations, other models, and reanalyses

Parameter sensitivity experiments show that the 4-yr mode of intrinsic variability is reasonably robust and appears whenever dissipation (both explicit and numerical) is sufficiently weak. For example, halving or doubling the GM coefficient does not significantly alter the properties of the Rossby waves in the Indian Ocean, although the amplitude of the waves is reduced at higher values of the GM coefficient. The Indian Ocean is subject to significant forced and coupled interannual variability, which could explain why this mode has not been previously noted. In particular, the period of the 4-yr mode is close to the 4–5-yr band of variability due to ENSO, so distinguishing between the two signals is likely to be difficult. Nevertheless, comparisons to observations and other ocean models are necessary to determine if this mode is due to a particularity of the model configuration or forcing.

Models forced by repeating annual forcing offer the best hope for clearly separating intrinsic from externally forced variability. Sérazin et al. (2015), using a suite of models, some with annually repeating forcing, found a band of enhanced, intrinsic, interannual variability in the south Indian Ocean stretching from Western Australia to Madagascar (the same region occupied by the 4-yr mode), but the detailed structure of and mechanism for this variability was not explored in their study. A similar band of enhanced variability in the south Indian Ocean was observed in a version of the ECCOv2 model (Wunsch et al. 2007) driven by annually repeating forcing (Piecuch and Ponte 2013); the propagation characteristics of this variability were consistent with mode-2 and mode-3 baroclinic Rossby waves (C. Piecuch 2016, personal communication).

AVISO (http://www.aviso.altimetry.fr) satellite altimetry shows variability in the south Indian Ocean with a period of 4 yr, but the SSH anomalies appear to propagate eastward and are thus unlikely to be associated with the westward-propagating 4-yr mode. It is perhaps not surprising that this intrinsic variability is difficult to identify in satellite altimetry; the satellite era is sufficiently short that only 4–5 oscillations could be present in the data, which is likely insufficient to separate the intrinsic variability from ENSO or make a good determination of phase.

The south Indian Ocean in the ECCOv4 state estimate (Forget et al. 2015), which uses the same grid as LLC–CNYF, hosts a westward-propagating disturbance with a 4-yr period that appears to originate from the west coast of Australia. However, this disturbance has a first-baroclinic-mode structure and thus propagates much more quickly than the intrinsic mode discussed in this paper. The LUC in ECCOv4 occupies a significantly larger depth range near the Australian coast than in LLC–CNYF, resulting in a substantially weaker core velocity and a more pronounced first-baroclinic-mode structure than in LLC–CNYF.^{6} Since the 4-yr mode identified in this paper appears to originate as an instability of the LC and LUC (see section 4), the vertical structure of this current system likely selects the vertical structure of the intrinsic mode. Observations of the LUC are sparse, so it is not clear which model more faithfully represents the LUC. The LUC transport in both ECCOv4 and LLC–CNYF is approximately 5 Sv—approximately the same as the observed value given by Thompson (1984)—but the current is horizontally larger than observed, resulting in peak velocities smaller than observed by more than a factor of 10 in both models.

The ECCOv4 product also suffers from the same short record length as AVISO. The estimation window of ECCOv4 is 20 yr, which is approximately the basin-crossing time of the intrinsic mode. If the intrinsic mode were not present in ECCOv4’s estimate of the initial conditions, it would not have time to fully develop within the estimation window, even if the model dynamics were conducive to the generation of the mode.

While it is hoped that an exhaustive search of observational products, ocean models, and reanalyses will lead to an unambiguous detection of the intrinsic mode discussed here, this undertaking is beyond the scope of the present paper and will be a subject of future research.

## 6. Summary

A coarse-resolution model forced by an annually repeating atmosphere supports an intrinsic mode of interannual variability that originates in the southern Indian Ocean and takes the form of multiple beams of long Rossby waves whose vertical structure is a mixture of the second and third baroclinic modes. This intrinsic variability also projects into the South Atlantic Ocean and so organizes variability across multiple ocean basins. The mode has a period of 4 yr with an amplitude and phase that vary on decadal time scales. Energy diagnostics suggest that this mode is driven by baroclinic instability of the Leeuwin Current and calculations using a 2.5-layer QG model show that this mechanism is plausible. The QG model also serves to illustrate that radiating baroclinic instability of an eastern boundary current generates waves with scales much larger than both the Rossby radius and the boundary current itself. This is contrary to the usual expectation that baroclinic instability produces disturbances with scales close to the Rossby radius.

While this 4-yr mode dominates interannual variability in the present model, similar forms of variability have only been identified in models where forced and coupled interannual variabilities have been suppressed by forcing the ocean with repeating atmospheric states. This fact may give the impression that the intrinsic mode discussed here is a mere curiosity. However, this analysis serves as a useful counter to the common perception that large-scale oceanic variability is primarily driven by interactions with the atmosphere. Unraveling the specific effects of such intrinsic variability on the predictability of the ocean and climate system is the subject of ongoing research.

## Acknowledgments

CLW, PC, and BDC are supported by the National Science Foundation under Grant OCE-1258887. Computational resources were provided by the Texas Advanced Computing Center and the San Diego Supercomputer Center through XSEDE. We also wish to thank two anonymous reviewers whose insightful comments improved this manuscript.

### APPENDIX A

#### Mechanical Energy Budget

Consider the inviscid, hydrostatic, Boussinesq primitive equations on a sphere, written in flux form

where the divergence is

**u** is the horizontal velocity vector, **v** is the 3D velocity vector, *ϑ* is the latitude, and *a* is Earth’s radius.

Let an overbar represent the time average and decompose the dynamical fields into mean and perturbation quantities, for example, . The disturbance flow evolves according to

and

Dotting the disturbance horizontal momentum *ρ*_{0}**u**′ into the horizontal momentum equations and averaging gives an evolution equation for the EKE :

where

is the instantaneous EKE, and the gradient in geographical spherical coordinates is

The unit vectors , , and point in the zonal, meridional, and vertical directions, respectively. Note that while the Coriolis terms cancel between the two momentum equations, the metric terms (proportional to ) remain since they are quadratic in the velocity.

The rhs of (A10) can be put in a more standard form by noting that

where the last equality follows from the continuity and hydrostatic equations. The continuity equation also allows the velocity correlation terms on the lhs of (A10) to be rewritten in terms of cross-gradient momentum fluxes:

Finally, the metric terms can be written in vectorial form using the identity

The final form of the EKE equation is

where the three source terms are

and

The term represents the conversion of potential energy into EKE, while the terms and represent conversion of mean kinetic energy into EKE by horizontally and vertically sheared flow, respectively. The term in proportional to arises from the metric terms in the horizontal momentum equations.

In studies using eddy-resolving models, the baroclinic conversion term is often further dissected using the buoyancy variance equation (e.g., Beckmann et al. 1994). Coarse-resolution models typically make use of the GM parameterization, which appears as an additional flux in the buoyancy equation and acts as a sink in the buoyancy variance equation. We avoid this complication by focusing on the kinetic energy budget, which does not involve GM. The velocities and eddy fluxes appearing in (3) and (4) are thus the velocities resolved by the model.

The conversion terms were calculated offline using snapshots of the velocity, temperature, and salinity. The shear production term was small compared to and and so is not shown. Note that a number of exact vector calculus identities are required to put the energy budget into standard form; the finite-difference schemes used in numerical models do not necessarily satisfy these identities exactly, although they are satisfied approximately. Thus, the conversion terms identified in (A18)–(A20) are likely not precisely equal to the conversion terms that would appear in an exact finite-difference EKE budget. Formulation of an exact EKE budget would be useful but beyond the scope of this study.

### APPENDIX B

#### Linear Stability Analysis

We examine the linear stability of a purely meridional mean current in the context of a 2.5-layer QG approximation. We consider a domain that is periodic in the *y* direction of length *L*_{y}. The longitudinal domain is −*L*_{x} ≤ *x* ≤ 0. The flow is governed by

where the potential vorticities in the two layers are

and

Here, *ψ*_{1,2} are the horizontal streamfunctions, and *β* is the meridional gradient of the Coriolis parameter. The squared wavenumbers *F*_{n}, defined as and , are related to the first and second Rossby deformation radii *L*_{1,2} through

where *L*_{1} and *L*_{2} correspond to the negative and positive signs, respectively, in (B4). The values of the layer depths *H*_{1,2}, the Coriolis parameter *f*, and the reduced gravities *g*′ and *g*″ are adjusted to approximately suit the conditions of the primitive equation computations (cf. Table B1). The meridional extent of the domain, which determines the meridional wavenumber, is chosen to maximize the growth rate of the unstable mode. We assume that the streamfunction is given by

with confined to the eastern boundary, located at *x* = 0. We take the form

The system (B1) is linearized around the mean meridional current, and the resulting equation for the perturbation is solved by assuming the form

The complex functions satisfy

subject to the boundary conditions

The radiation condition (B12) is determined by matching the solution of (B8)–(B11) for *x* > −*L*_{x} to free Rossby waves [i.e., the solution to (B8)–(B11) with *V*_{1} = *V*_{2} = 0 for *x* < −*L*_{x}] with a given frequency *ω* and meridional wavenumber *l*. The free waves must have westward group velocity and remain finite as *x* → −∞, which requires the imaginary part of the wavenumber to satisfy *k*_{i} ≤ 0. For growing or neutral waves, these criteria select free long Rossby waves that have the wavenumbers

for the first and second baroclinic modes, respectively. The coefficients and are then the projections of *ϕ*_{1,2} at *x* = −*L*_{x} onto the two free Rossby waves having wavenumbers *k*^{(1)} and *k*^{(2)} (cf. Hristova et al. 2008).

The system (B8)–(B12) amounts to a set of homogeneous, coupled, ordinary differential equations where the complex frequency *ω* is the eigenvalue and *ϕ*_{1,2}(*x*) are the complex eigenfunctions. This system is solved numerically, with a grid that resolves the smallest of the two deformation radii. An example of a typical unstable solution, with a meridional wavenumber *l* chosen to be the gravest mode fitting the domain, is shown in Fig. 14.

## REFERENCES

*Fourier Analysis of Time Series: An Introduction.*2nd ed. Wiley, 288 pp.

*Descriptive Physical Oceanography: An Introduction.*6th ed. Academic Press, 555 pp.

*Statistical Analysis in Climate Research.*Cambridge University Press, 484 pp.

*Ocean Circulation and Climate*, G. Siedler, J. Gould, and J. Church, Eds., Academic Press, 47–58.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JPO-D-16-0177.s1.

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).

^{1}

Several values between 250 and 1250 m^{2} s^{−1} were also tried, but their climatologies were less satisfactory.

^{2}

This would be at a confidence of 95% relative to white noise.

^{3}

Estimated using the WKB approximation (see, e.g., Chelton et al. 1998) and averaged over the box.

^{4}

Note that the phase is approximately −30° for the last 60 yr of the time series, and it is likely that regular phase shifts of ~90° would continue to occur if the time series were extended past 300 yr. The approximate bimodality of the time series is thus probably an artifact of its length, and the true phase distribution is multimodal.

^{5}

The projection was determined by regression of the real and imaginary parts of CEOF 1 on *δ***x**.

^{6}

This may be due to the larger value of the GM parameter and the more diffusive advection scheme used in ECCOv4 as well as to differences in Southern Ocean buoyancy forcing which impact the formation of the water masses that make up the Indian Ocean.