## Abstract

Low-frequency variability in the south Indian Ocean is studied by analyzing 200 years of output from a fully coupled climate model simulation. At time scales of 2–10 years, the variability is dominated by westward-propagating features that form on the eastern side of the basin. Using feature tracking and clustering, the spatiotemporal characteristics and preferred pathways of the propagating features are identified and studied in detail. By comparison of the phase speed and vertical structure of the propagating anomalies identified by the feature tracking with linear theory, we conclude that these features are likely mode 1 or 2 baroclinic planetary waves. The effects of this low-frequency variability on the climate system is investigated. By analysis of the mixed-layer temperature budget, it is shown that at particular geographic locations, the propagating features can substantially modify the near-surface ocean and induce significant fluxes of heat into the atmosphere. In turn, these heat fluxes can drive a coherent atmospheric response, although this response does not appear to feed back onto the ocean. Finally, we discuss the implications for the interannual climate predictability.

## 1. Introduction

Understanding the low-frequency variability of the coupled ocean–atmosphere system on time scales of 2–10 years, commonly referred to as “near term” variability, presents both a major challenge and an exciting opportunity. The challenges arise due the fact that variability within the climate system is a complex, multiscale problem: “fast” processes, such as atmospheric forcing or oceanic mesoscale turbulence, can excite or modify slower processes such as large-scale variations in the thermocline (Hasselmann 1976; Franzke et al. 2015). As such, a dynamical understanding of low-frequency variability in the ocean–atmosphere system requires an understanding of the complex interactions between phenomena with a large range of temporal and spatial scales.

The ocean is capable of supporting propagating disturbances and teleconnections, which provides a useful framework for understanding low-frequency variability within the broader climate system (O’Kane et al. 2014). Since the time scales are considerably longer in the ocean than in the atmosphere, oceanic teleconnections can provide a mechanism for interseasonal and interannual predictions. Of particular interest are long baroclinic planetary waves (LPWs), which have the following characteristics (Killworth et al. 1997; Killworth and Blundell 1999; Maharaj et al. 2007):

wavelengths that are long relative to the local deformation radius

*R*_{D};periods that are long relative to the inertial period

*f*^{−1}of ~24 h;scale separation, both spatially and temporally, between the perturbations due to planetary waves and the climatological state.

LPWs propagate slowly and remain relatively coherent over long time scales (White 2000; Fu and Qiu 2002). As such, they are able to transmit energy and information over large distances. They have clear signatures in sea surface height (Chelton and Schlax 1996; Morrow and Birol 1998; Birol and Morrow 1999; Fu and Qiu 2002; Fu 2004; Maharaj et al. 2007), salinity (Menezes et al. 2014; O’Kane et al. 2016), ocean color (Cipollini et al. 2001), and, of particular relevance to this study, sea surface temperature (Halliwell et al. 1991a; Cipollini et al. 1997; White et al. 1998; White 2000; White et al. 2004). LPWs manifest at depth as significant undulations of the thermocline, which have been detected in situ from long-term XBT data (White and Saur 1981; Kessler 1990; Masumoto and Meyers 1998) and ocean reanalyses (White et al. 2004; Vargas-Hernandez et al. 2015). They are found in all major ocean basins and provide a route by which information may propagate from the midlatitudes into the tropics (O’Kane et al. 2014; Wolfe et al. 2017).

Despite the fact that the basic characteristics of LPWs have been extensively studied, their role within the climate system and influence on near-term climate predictions are relatively unknown. There is some evidence from both observations and numerical models that propagating LPWs can induce persistent sea surface temperature anomalies. In particular, the work of Frankignoul (1981) using mooring data, and of Halliwell et al. (1991a,b) and Leeuwenburgh and Stammer (2001) using satellite-derived sea surface temperatures (SSTs) and sea level anomalies, showed that LPWs were able to redden the SST signal and enhance its variance compared with SST forced by fast atmospheric motions alone. Theoretical studies employing idealized physics have shown that slow LPWs can have a substantial impact on SSTs by driving entrainment of cooler thermocline water into the mixed layer or through anomalous advective heat fluxes (Stevenson 1983; McGillicuddy and Robinson 1998; Goodman and Marshall 1999). Simple coupled models have also demonstrated that the SST anomalies induced by LPWs can drive anomalous winds that can feed back onto SST anomalies (Hasselmann 1976; Frankignoul and Hasselmann 1977; White et al. 1998). A positive feedback can occur in certain parameter regimes, where the local atmospheric response acts to amplify the SST perturbations which, in turn, deepen the subsurface thermocline perturbations (White et al. 1998; Goodman and Marshall 1999; White 2000).

In this study, we use a 200-yr-long run of a fully coupled climate model to study oceanic teleconnections due to LPWs and their direct influence on the sea surface temperature and the near-surface atmosphere. We focus on the subtropical south Indian Ocean because of its role in the global climate and also because it is comparatively understudied when compared with the North Atlantic or the Pacific. However, the state of the south Indian Ocean can have far reaching consequences, influencing, for example, rainfall variability in East Africa, Madagascar, and Indonesia (Schott et al. 2009) and the primary productivity and fish abundance off the African coast (Xie et al. 2002; White et al. 2004). As such, there is a strong societal demand for improved near-term prediction capability, which requires a better understanding of the dynamics and variability in this region. That said, the approach taken is quite general and could be easily applied to other geographic regions.

Previous work investigating variability in the south Indian Ocean has revealed that the region is host to variability on a variety of time scales (Schott et al. 2009; Han et al. 2014) and that LPWs strongly influence the near-surface ocean on interannual time scales (Masumoto and Meyers 1998; Birol and Morrow 1999; White 2000; White et al. 2004). In particular White (2000), using 6 years of satellite sea surface height (SSH) and SST data, complemented by 29 years of atmospheric and SST reanalysis, found clear evidence of slowly propagating subtropical LPWs with periods of 2–5 years, emanating from the coast of Australia and propagating westward at phase speeds of between 4 and 7 cm s^{−1}. He also found in phase covariation of the SST and SSH anomalies, and evidence of coupling between the oceanic planetary waves and atmospheric winds, although the strength of this coupling varied strongly with time. By careful examination of the phase relationships between the SSH, SST, and surface wind anomalies, White (2000) concluded that the loss of heat from the upper ocean to the atmosphere induced by the LPWs was approximately balanced by the anomalous heat gain by the suppression of the entrainment of cooler thermocline waters into the surface mixed layer. In a follow-up study, White et al. (2004) noted covariance between the surface fields and deeper, pycnocline anomalies, and that the passage of LPWs might influence the fish stocks of the region.

However, recent work on south Indian Ocean LPWs has revealed complicated physical and propagation characteristics. Using a coarse-resolution forced ocean general circulation model, O’Kane et al. (2014) revealed that midlatitude LPWs propagate along a limited set of favored pathways. In the south Indian Ocean, these pathways were found to be orientated nonzonally, and the features propagating along them underwent substantial changes in shape, amplitude, and phase speed as they traveled, including strong interactions with bathymetry. In a follow-up study (O’Kane et al. 2016), it was shown that this mode of variability was present in the surface salinity fields of a model forced with both repeat year and interannually varying atmospheric forcing. In a similar study, Wolfe et al. (2017) investigated LPWs originating from the west coast of Australia in a model forced by repeat year atmospheric forcing. The LPWs isolated by Wolfe et al. (2017) propagate across the Indian Ocean Basin along multiple nonzonal pathways, before either impacting the African continent or leaking into the South Atlantic. When LPWs reach western boundaries, they can induce large modifications in the transport, stratification, and eddy structure of the associated boundary currents and in the Agulhas Current system (Schouten et al. 2002a,b; Wolfe et al. 2017). The fact that both O’Kane et al. (2016) and Wolfe et al. (2017) found similar phenomena in models forced by repeat year atmospheric forcing suggests that this mode of variability is intrinsic to the Indian Ocean and not necessarily excited by low-frequency variability in the forcing itself.

In this paper, we demonstrate that our numerical model produces a mode of low-frequency variability that manifests as slowly propagating anomalies, reminiscent of the baroclinic LPWs found by White (2000). These anomalies follow distinct pathways defined by strong gradients in potential density, as described in O’Kane et al. (2014). We are able to isolate these distinct pathways by combining feature tracking with clustering analysis, and to characterize how the nature of the propagating anomalies differs between the different pathways. We apply a mixed-layer temperature budget to quantify the response of SST to the low-frequency oceanic variability, and the potential impact of the anomalies on the near-surface atmospheric circulations. In a departure from previous work on the subject, we do not find clear evidence for strong atmospheric feedback. Finally, we will discuss the implications of work for predictability on interannual time scales. Our work builds and expands upon the results of O’Kane et al. (2014) by combining their approach with elements of the studies of Halliwell et al. (1991a,b), Leeuwenburgh and Stammer (2001), and White (2000).

The remainder of this article is organized as follows. Section 2 briefly introduces the background theory of LPWs and the framework used to quantify their impact on sea surface temperatures. Section 3 describes the numerical model configuration used in this study. Section 4 describes in detail the nature of the low-frequency variability in the Indian Ocean Basin. The influence of the propagating anomalies on the ocean state, and particularly the SST, is discussed in section 6. Finally, in sections 7 and 8 we discuss the implications of our work, with a particular focus on potential for the ocean variability to drive a coherent atmospheric response, and the implications of this study for near-term climate predictions.

## 2. Theoretical background

### a. Free, forced, and coupled long planetary waves

Classical Rossby wave theory is developed in numerous textbooks (e.g., Gill 1982, 159–175). Here, we present a very brief overview to orient the reader and to provide a framework for interpreting our results. A reader seeking more details is referred to the textbook cited above. Consider a quasigeostrophic, continuously stratified, flat-bottomed ocean with a background buoyancy-frequency profile *N*^{2}(*z*) = −(*g*/*ρ*_{0})*dρ*/*dz*. By linearizing the momentum and continuity equations, projecting onto vertical normal modes, and using a plane wave trial solution, the homogeneous (unforced) problem yields the Rossby wave dispersion relationship for the *m*th normal mode:

where *ω* is the temporal frequency, *k* and *l* are respectively the zonal and meridional wavenumbers, and *R*_{m} is the *deformation radius* for the *m*th mode. The normal modes, *ϕ*^{m}(*z*), and the deformation radius, can be obtained by the solution of the linear Sturm–Liouville eigenvalue problem (Smith 2007):

The mode for which *m* = 0, called the barotropic mode, gives rapidly propagating Rossby waves and is not of interest in this study. The higher modes *m* > 0, known as baroclinic modes, propagate substantially more slowly than the barotropic mode.

In the long-wave limit of $k,l\u226aRm\u22122$, often referred to as the planetary wave limit, the phase speed of the *m*th mode is

Hence, the phase propagation of linear LPWs in the absence of topography or background flow is nondispersive and set entirely by the background stratification profile *N*^{2}(*z*) and *β*. Additionally, we note that Eq. (3) is everywhere less than 0, implying *westward* phase propagation. An extended theory of extratropical LPWs described in Killworth et al. (1997) and Killworth and Blundell (1999, 2003) has shown that both bottom bathymetry and the mean flow can also influence the speed and direction of propagation, although this extended theory does not completely explain the observed dispersion characteristics (Maharaj et al. 2007).

The dispersion relationship in Eq. (1) is derived for *free* Rossby waves, as the role of forcing has been excluded by treating only the homogeneous problem. In the real ocean, however, free Rossby waves exist concurrently with *forced* Rossby waves that arise due to local variations in the wind stress curl (White 1977; Masumoto and Meyers 1998). As such, variability at any given location is likely to arise from a mixture of local forcing and the influence of waves generated remotely. White (1977) has shown that, in cases of resonance between the free waves and local forcing, the phase speed of the forced waves is twice that of the free wave given by Eq. (3). More detailed analysis by Qiu et al. (1997) suggests that although forced Rossby waves generally move at phase speeds greater than the free waves, blending of free and forced waves can result in behavior that reinforces the apparent propagation at the long-wave phase speed. As such, disentangling observed free and forced solutions is difficult and will not be attempted here. However, in general the presence of forcing results in waves that propagate faster than the theory of free waves would suggest. This point will be of relevance when we consider the phase speed of propagating features in our mode.

Finally, we note that LPWs can couple with the atmosphere to create a self-sustaining feedback whereby the SST variations induced by planetary waves can then modify the heat flux between the atmosphere and ocean, and hence the local atmospheric circulation. In certain parameter regimes, the anomalous ocean circulation can exert an anomalous stress that amplifies oceanic Rossby waves (Goodman and Marshall 1999). White (2000) and White et al. (2004) found evidence of these feedback processes in the south Indian Ocean. By employing a simple mixed-layer temperature budget (see section 2b), and assuming a linear relationship between the meridional wind anomaly and the SST anomaly, White (2000) was able to derive a modified dispersion relationship that incorporated the effects of this coupling. Numerical estimates of the coupled phase speed in the Indian Ocean indicate that the coupled Rossby waves propagate westward at between 1/2 and 1/3 of the baroclinic mode 1 free Rossby wave speed. This is in contrast to forced Rossby waves, which generally travel westward faster than free waves.

### b. The ocean mixed-layer temperature budget and the low-frequency variability of SST

In this study, the basic approach used to separate the influence of slow oceanic motions from that of high-frequency atmospheric forcing is a simplified, one-dimensional mixed-layer temperature budget. Here, we present a sketch of the derivation. The full derivation is included in appendix B. We begin with the equation for the conservation of conservative temperature Θ, written in flux form:

where ∇_{xy} = (∂/∂*x*, ∂/∂*y*) is the horizontal gradient operator, **u** = (*u*, *υ*) is the horizontal velocity, *w* is the vertical velocity, *q* is the heat flux, and the rest of the notation is standard.

We now assume a homogeneous mixed layer with a time-varying depth *h*. Integrating from the surface at *z* = 0 (ignoring variations of the free-surface) to *z* = −*h* and assuming that all dependent variables with the exception of the vertical velocity are constant within the mixed layer, such that for some quantity *φ*, then $\u222b\u2212h0\phi dz=h\phi $, we obtain

where ΔΘ = Θ − Θ^{+} is the temperature difference between the mixed-layer temperature and the temperature just below the base of the mixed-layer Θ^{+}, and *w*_{e} is the *entrainment velocity*, given by the kinematic condition at the base of the mixed layer:

Note that *w*_{e} can only be positive or zero, as the loss of fluid from the mixed layer into the interior (detrainment) cannot change the mixed-layer temperature. Also, *Q* is the total surface heat flux (the sum of shortwave, longwave, sensible, and latent heating of the ocean surface). In deriving Eq. (5) we have ignored heat fluxes below the mixed layer. We then perform a Reynolds decomposition, separating Θ, **u**, and *h* into a seasonally varying mean $(\u22c5)\xaf$ and fluctuating component (·)′. As the mixed-layer temperature budget consists of triple product terms, the resulting decomposition is complicated. Further details of the derivation are presented in appendix B. The Reynolds decomposition allows us to derive an expression for the mixed-layer temperature anomalies:

Note that the advective and entrainment terms, as well as the product of the mixed-layer depth anomalies and the temperature tendencies, operate on the full temperature field (i.e., $\Theta =\Theta \xaf+\Theta \u2032$). Following Frankignoul (1985), we note that the mean entrainment term $w\xafe\Delta \Theta \u2032$ is normally positive during the cool season when the mixed layer deepens, and zero otherwise. Assuming that the anomalous thermocline temperature is small $(\Theta \xaf+\u2032\u22480)$ on seasonal time scales, then the mean entrainment acts to damp SST anomalies exponential with inverse decay time scale of $we\xaf/h\xaf$ [see Eq. (15) in Frankignoul (1985)]. As such, we write

where *λ* is the (inverse) damping time scale. We now interpret the various terms in Eq. (7) by first decomposing the horizontal velocity into a geostrophic component, denoted by the subscript *g*, and an ageostrophic component (subscript ag) the latter of which is dominated by the wind-induced Ekman currents such that **u**_{ag} ≈ **u**_{Ek} given by

where *τ*_{wind} is the surface mechanical wind stress. Equation (8) relies on the assumption that the Ekman layer depth is less than *h*, which is justifiable in the subtropics and midlatitudes, but may not be applicable in the tropics. As such, the contribution of Ekman divergence may be overestimated in tropical regions. Equation (7) can now be written as follows:

where

is the forcing due to the surface heat flux,

is the wind forcing, and

is the forcing due to ocean processes processes.

The nonlinear mixed-layer temperature budget in Eq. (9) allows us to diagnose the relative impact of anomalous surface heat fluxes, wind forcing, and purely oceanic effects such as anomalous temperature advection and entrainment on the near-surface temperatures. However, we note that the mixed-layer temperature budget analysis is unable to explicitly separate the influence of LPWs from other dynamical processes. Instead, we infer the influence of LPWs by a detailed analysis of the time scales present in the various terms of the mixed-layer temperature budget and the SST. To illustrate, we take the Fourier transform of Eq. (9). Assuming that the mean advection term is small, we obtain

where $(\u22c5)\u02dc$ denotes the Fourier transform. Noting that the power spectral density (PSD) *P*_{xx} of a variable *x* is given by $Pxx(\omega )=x\u02dcx\u02dc*$, we square each side of Eq. (13) to give an expression for the temperature PSD:

where $Cxy(\omega )=X\u02dc(\omega )Y\u02dc*(\omega )$ are the cross-spectra for two variables *x* and *y*, and the operator Re indicates that the real part of the complex cross spectra should be taken. In Eq. (14) the subscript *W* represents the wind forcing, *O* the oceanic forcing, and (as previously defined) *Q* the total surface heat flux.

We now consider the PSDs and cross spectra of the forcing terms of the right-hand side of Eq. (14) and to be stochastic inputs. If there is no anomalous oceanic forcing, and if we take the wind forcing and surface heat flux to be a process with a constant power spectrum [i.e., white noise, as described in Frankignoul (1985), Hall and Manabe (1997), and Leeuwenburgh and Stammer (2001)] then Eq. (14) reduces to an order-one autoregressive process [AR(1)] with the spectra

where *A* = constant. At low frequencies (i.e., where *ω*^{2} ≪ *λ*) the temperature power spectrum is flat, whereas at high frequencies (where *ω*^{2} ≫ *λ*) the temperature power spectra resemble a red-noise process with a spectral roll off *ω*^{2}. The value of the damping coefficient *λ* controls the transition between these regimes (Nise 2011, 545–547).

However, in the presence of oceanic forcing with periods longer than the inverse damping time scale and (for illustrative purposes) negligible cross-spectra between the atmospheric and oceanic forcing, the SST spectrum resembles the power spectra of the slow oceanic process (scaled by the inverse square of the damping coefficient) at low frequencies, that is,

As such, if the oceanic forcing given by Eq. (12) has an influence on the mixed-layer temperature at low frequencies, we expect both the oceanic forcing and the temperature power spectra to have peaks in similar frequency bands. Equivalently, if the influence of oceanic forcing is weak, we expect a flat power spectrum at low frequencies.

Applications of this approach using both observations and numerical models, similar to those presented in this study, can be found in Frankignoul (1981), Frankignoul and Reynolds (1983), Halliwell et al. (1991a,b), and Leeuwenburgh and Stammer (2001).

## 3. The coupled climate model

In this paper, we use a 500-yr-long run variant of the GFDL Climate Model 2.1 (CM2.1) (Delworth et al. 2006), used within the CSIRO Climate Analysis Forecast Ensemble (CAFE) near-term prediction system (O’Kane et al. 2019). The model uses the same atmospheric, land, and sea ice components as CM2.1 (AM2, LM2, and SIS, respectively). However, the ocean uses a newer ocean model (MOM4p1). The ocean model grid is the tripolar ACCESS-o grid (Bi et al. 2013) with a nominal grid spacing of 1° but with a finer latitudinal grid spacing in the tropics and the Southern Hemisphere high latitudes, (telescoping to 0.33° at the equator and to 0.25° in the Southern Ocean at 75°S). There are 50 vertical levels, with 10-m grid spacing in the upper ocean, increasing to a maximum of 300 m. Subgrid processes for the ocean model are adopted from CM2.1, including neutral physics (Redi diffusivity and Gent–McWilliams skew diffusion), Brian–Lewis vertical mixing profile, Laplacian friction scheme, and a K-profile parameterization for the mixed-layer calculation. The atmospheric model (AM2) has a grid spacing of 2° latitude and 2.5° longitude, and 24 hybrid (sigma–pressure or terrain following pressure) vertical levels. Concentrations of atmospheric aerosols and radiative gases, and land cover are based on 1990 conditions.

To reduce model biases in the properties of mode waters, the model temperature and salinity fields are restored to *World Ocean Atlas 2013* (*WOA13*) climatology at depths below 2000 m, with a restoring time scale of 1 year. This deep restoring improves the model’s representation of the upper ocean stratification, at the expense of suppressing variability with multidecadal time scales. However, we note that this study is concerned with variability on much shorter time scales (from 1 to 10 years) and that variability on longer time scales is not considered in this paper. We have tested this assumption by applying some of the diagnostics used in this paper (notably, the in-band variance described in section 4a) to a similar model with no deep restoring. Minimal differences in the upper ocean variability in Indian Ocean are found. An additional benefit of the deep restoring is that the spinup time is significantly reduced when compared with the free model, which after 500 years of spinup is still drifting significantly.

The model is run for 500 years, with the first 300 years used to spin the model up from initial conditions set by the *WOA13* climatology. Only the final 200 years of the simulation are used for analysis. The oceanic fields used in this study are saved at monthly intervals, while atmospheric fields are saved daily.

## 4. Low-frequency variability in the Indian Ocean

### a. Characterization of the interannual mode of variability

#### Variance in discreet time bands

To characterize the spatial patterns of variability within the Indian Ocean Basin, we calculate the variance of conservative temperature within specified temporal bands (henceforth referred to as the *in-band variance* (Monselesan et al. 2015; O’Kane et al. 2016). Specifically, to determine the variance contained in all spectral bands of a desired variable between frequencies *f*_{1} (with corresponding period *T*_{1}) and *f*_{2} (period *T*_{2}), at each model grid point, we compute the power spectral density (PSD) from the final 200 years of our simulation using Welch’s method, then integrate the resulting PSD from *f*_{1} to *f*_{2} numerically using the trapezoidal rule that, by Parseval’s theorem, gives the resulting variance.

The in-band variance of the conservative temperature anomaly (defined relative to a seasonally varying mean) is computed for frequency bands corresponding to periods of less than 1 year, 1–2 years, 2–5 years, and 5–10 years and shown in Fig. 1. At the surface the temperature variance is pronounced north of 10°S at higher frequencies (periods of less than 2 years) (Figs. 1a,d), with lower-frequency variance concentrated between about 30° and 10°S (Figs. 1g,j). At 150-m depth (Figs. 1b,e,h,k), the high-variance regions are organized along several beams. At high frequencies (period < 2 years), high variance is found at approximately 10°S, running zonally from south of the island of Java to the northern tip of Madagascar, following the path of the South Equatorial Current. At lower frequencies (period > 2 years), variance remains high in the beam located at 10°S, indicative of a broadband signal with power distributed across a range of different frequencies. However, there is also significant variability in the subtropics: a nonzonal region of high variance extends northwestwards from the southern tip of Western Australia to approximately 20°S, 60°E. A weaker, more zonally oriented region is found at approximately 30°S, extending from the coast of Australia to 80°E. In all cases, these beams broadly align with regions of high potential density gradients (solid lines in Fig. 1), suggestive of waveguides.

A latitude–depth section at 90°E (Figs. 1c,f,i,l) shows the vertical structure of this mode of variability. At this longitude, the variance at time scales longer than 2 years is found in two distinct regions, with maxima found at approximately 27° and 15°S, and at depths of between 100 and 150 m, although elevated variance is found as deep as 450 m. In contrast, in the higher-frequency temporal bands (less than 2 years) the high variance is closer to the surface for the northern beam, and is absent below about 200 m.

### b. Pathways of propagating anomalies

In the previous section we have established that variability on 2–10-yr time scales is organized into a series of narrow regions that cross the Indian Ocean. Complex empirical orthogonal function analysis undertaken by O’Kane et al. (2014) and Wolfe et al. (2017) provides evidence that the high variance in these regions is due to the presence of propagating features that follow a set of preferred pathways. However, both O’Kane et al. (2014) and Wolfe et al. (2017) have shown that pathways are both nonzonal and depth dependent, results that are confirmed by our in-band variance calculations. This poses a challenge: How can we best define propagation pathways that evolve in three dimensions?

To attack this problem, we first remove the depth dependence by computing the geostrophic dynamic height *ψ*_{g}, given by

where *S*_{A} is the absolute salinity, *p* is the pressure, *S*_{O} = 35.165 04 g kg^{−1} is the reference absolute salinity, $\nu ^$ is the specific volume, and *p*_{ref} is the reference-level pressure (McDougall and Klocker 2010). The dynamic height anomaly is defined such that

that is, the horizontal gradient of *ψ*_{g} gives the geostrophic current velocity relative to the current at a reference pressure *p*_{ref}. The advantage of this approach is that *ψ*_{g} is a vertically integrated quantity. As such, the value of *ψ*_{g} at the surface is a function of all values of temperature and salinity between the surface and the reference depth. By using the dynamic height, we transform a three-dimensional problem into a two-dimensional problem.

We compute *ψ*_{g} at all model grid points and time steps in the 200-yr simulation, using the TEOS-10 Gibbs Seawater software package (McDougall and Barker 2011). We define the anomalous *ψ*_{g} relative to a monthly climatology with *p*_{ref} set to 500 db. The 500-db level is chosen as it is just below the depth of the high variance core in the interannual time bands, as indicated by Figs. 1c–f, filtering the influence of variability deeper than the features of interest. As a result, the propagating features studied here appear more clearly in the 500-db referenced dynamic height field than in the sea surface height field (which is equivalent to setting *p*_{ref} to the ocean depth) determined directly by the model.

To objectively define the propagation pathways, we apply a feature detection and tracking algorithm to the surface dynamic height anomalies. Here, we adapt the method described in Chelton et al. (2011) for the detection of mesoscale eddies from sea surface height altimetry. Although the analogy between turbulent mesoscale eddies and LPWs may not be immediately apparent, we justify the use of this methodology by first noting that both eddies and LPWs manifest as contiguous regions of anomalous dynamic height enclosed by a single contour. The feature detection algorithm first identifies regions of dynamic height anomaly of one or more extrema (maxima in the case of anticyclones, minima in the case of cyclones) that can be enclosed within a closed contour. The feature tracking algorithm then tracks these detected features through time by matching the individual features identified. The methodology is described in further detail in appendix A, together with a detailed discussion of individual steps of the algorithm and the justification for our choice of hyperparameters. The implementation is taken from the open-source software accompanying Oliver et al. (2015).

The preferred propagation pathways of these large-scale features are determined by a clustering analysis of the latitudes and longitudes of all detected features with lifetimes greater than 1 year. The clustering methodology groups the trajectories of these long-lived features into a number of distinct paths. Here, we use a Gaussian mixture model (GMM) to carry out the clustering (Press et al. 2007, 842–850). Unlike other clustering methods more common in climate science, such as K-means, GMMs permits clusters with anisotropic distributions in the original data space, which allows for the definition of clusters with a preferred orientation in geographic space. By exploiting the anisotropic nature of detected trajectories, we are able to define zonally elongated clusters. To determine the optimal number of clusters, we performed the clustering using 1 to 20 different clusters. We then used the Akaike information criterion (AIC), a statistical metric that balances the capacity of a (statistical) model to fit the data against the number of free parameters to guard against overfitting, to determine an “optimal” number of clusters (Burnham and Anderson 2004).

The clustering analysis, together with the AIC assessment, places the trajectories into five distinct clusters, shown in Figs. 2a and 2b. The five clusters (from north to south) are as follows:

One cluster (colored cyan) is oriented zonally at latitude 10°S and extends across the Indian Ocean Basin, terminating near 60°E. Features assigned to this cluster propagate westward and originate either in the shallow seas north of the Australian continent, or from disturbances formed south of the Indonesian archipelago.

Two clusters (colored green and red) originate off the Western Australian coast (near 20° to 35°S) and extend across the Indian Ocean Basin. Both these “plumes” have a substantial nonzonal orientation. Features in both these clusters propagate westward.

One cluster (colored purple) in the center of the basin at latitude of ~28°S is zonally oriented and arises from features that originate near longitude 90°E. Features assigned to this cluster propagate westward, dissipating when they impact the southern tip of Madagascar.

One cluster to the south of the domain (colored blue) originates near the coast of Africa and extends zonally about 25° of longitude. Features detected in this cluster propagate eastward.

To determine the mean pathways, we calculate the latitude centroid at each longitude from all points in each cluster. We then smooth the resulting path with using a 5°-wide Hamming filter (black dots in Figs. 2a and 2b). The mean pathways pass through the approximate center of the cloud of points associated with each cluster, and align roughly with the locations where features are more frequently detected.

In addition to determining the propagation pathways, the feature tracking also allows us to calculate the propagation (phase) speed by following the individual tracked features through time. In Fig. 2c we show both the mean (solid circles) and the two standard deviation limit (bars) in 1° latitude bins. Phase speeds are typically less than 10 cm s^{−1}, even in the tropical region, with an expected poleward reduction in phase speed. However, there is substantial variability among the speed of the detected features, with the 2 standard deviation level ranging from around 3 cm s^{−1} in the south to as much as 15 cm s^{−1} in the tropics.

Together with the phase speeds calculated from the feature tracking analysis, we also plot in Fig. 2 the theoretical LPW phase speed from the long-wave approximation in Eq. (3) for the first three baroclinic modes. To determine the deformation radius *R*_{m}, we solve the linear eigenvalue problem [Eq. (2)] numerically using the finite difference discretization of Smith (2007) at each model grid point. The background stratification in Eq. (2) is taken to be the time-mean in situ density from the final 200 years of the model run. Remarkable agreement is found between the empirically derived phase speeds and the mode 2 theoretical phase speed. In contrast, the first baroclinic mode is far too fast to describe the data, as it either follows or, at times, exceeds the 2 standard deviation envelope of the measured phase speeds.

The apparent similarity between the observed phase speed of the tracked features and the estimated theoretical second baroclinic mode phase speed requires further investigation. To this end, we follow O’Kane et al. (2014) and Wolfe et al. (2017) and calculate the complex empirical orthogonal functions (CEOFs), described in von Storch and Zwiers (1999, 353–365). In contrast to standard EOFs, CEOFs are capable of revealing the structure of *propagating* signals, provided that they are nondispersive. A CEOF has both a real and imaginary components, both in their spatial patterns and their associated principal component (PC) time series. Propagating components are indicated by real and imaginary components that are in approximate phase quadrature—that is, out of phase by 90°.

We compute the CEOFs using the three-dimensional conservative temperature anomaly field in the Indian Ocean Basin, weighted by the volume of each grid box. The first three CEOFs are shown in left, center, and right columns of Figs. 3; together they explain ~35% of the total variance. The eigenvalue spectrum is quite flat, with no clear separation between dominant and subdominant modes. However, the phase quadrature characteristic of propagating modes becomes less apparent in the PC time series for modes greater than three. As such, we limit our discussion here to the first three CEOF modes.

The horizontal spatial patterns of variability as determined by the first three CEOFs are shown in Figs. 3a(i)–(iii). There are two dominant spatial patterns evident in all three modes: north of 20°S, there is a multipole structure in both real and imaginary CEOFs. The first CEOF [Fig. 3a(i)] shows a bipolar modal structure with negative values lying to the south of Sumatra and positive values lying through an elongated region that traverses the basin zonally. The multipolar structure of variability north of 20°S aligns with the northern “beam” of high variance in Fig. 1. South of 20° there is evidence of a series of positive and negative anomalies originating near the coast of Western Australia. These patterns of variability are broadly replicated by the second [Fig. 3a(ii)] and third [Fig. 3a(iii)] CEOF modes. The PC time series, shown in Fig. 3, also show approximate quadrature between real and imaginary components, indicative of propagation, as well as quasi-periodic behavior with a time scale of 4–8 years.

A longitude–depth slice of the first three CEOFs along the southern propagation pathway, shown as a solid line in Figs. 3b(i)–(iii), gives some insight into the vertical structure of these propagating features. Much like in the plan view presented in Fig. 3a, real and imaginary components of the CEOFs are generally out of phase by between 90° and 180° longitude, and as such are in approximate quadrature. However, there is no strong evidence of a phase shift between real and imaginary components in depth, indicating that any vertical propagation of the temperature anomalies through the water column is weak. In all three CEOFs, the disturbances are confided to the upper 1000 m of the ocean. For a given phase, both the real and imaginary CEOFs have 1 to 2 zero crossings, reflective of low-order baroclinic modes. Indeed, projection of the CEOFs onto baroclinic modes determined by the numerical solution of Eq. (2) suggests that the first CEOF is dominated by baroclinic mode 1, whereas CEOFs 2 and 3 are dominated by mode 2.

The CEOF analysis provides further evidence that the propagating features detected by the feature tracking methodology can be considered to be LPWs, with a vertical structure and propagation characteristics dominated by baroclinic mode 2. This result is in contrast to the results of White (2000) and White et al. (2004), who suggested that the propagating modes in the Indian Ocean were likely baroclinic mode 1 features whose propagation speed was retarded by coupling with the atmosphere. Instead, our results support those of Wolfe et al. (2017), who attributed these propagating features to intrinsic ocean processes dominated by baroclinic modes 2 and 3. Additionally, the slow propagation of these waves strongly suggests that local forcing plays a minor role, as theoretical calculations suggest such waves should propagate more quickly than free waves (White 1977; Qiu et al. 1997).

Of the five clusters described above, only the first three will be considered further. The latter two clusters are associated with Subtropical Mode Water production (purple) and the Agulhas retroflection (blue). As the dynamics of features assigned to these clusters are distinct from those found in the three more northern clusters and the pathways do not connect the subtropics to the tropics, they are not considered further in this paper.

## 5. Characteristics of propagating anomalies

Having determined that the variability within the south Indian Ocean on 2–10-yr time scales is dominated by westward-propagating anomalies generated in the east, we now investigate the spatiotemporal characteristics of the anomalies as they move along these pathways.

Forty-year-long time–longitude Hovmöller plots of the dynamic height anomaly *ψ*_{g} at the ocean surface [computed from Eq. (17)] for the southern (Fig. 4a, top), central (Fig. 4b, top), and northern (Fig. 4c, top) pathways determined from the feature tracking calculation described in section 4b show clear evidence of propagation along all three pathways. The amplitudes of the anomalies are also roughly the same for each pathway (about 20 cm). However, there are differences in the propagation speed of the anomalies among the pathways. On the southern pathway (Fig. 4a, top) the anomalies propagate at varying speeds, taking an average of 6 years to move from the West Australian coast (110°E) to approximately 85°E. West of 85°E, the anomalies then accelerate rapidly, traversing the remaining 25° of longitude in less than 1 year. The location where this acceleration occurs is also coincident with a significant topographic feature: the Ninety-East Ridge. Similar behavior in the Pacific Ocean was attributed to resonant interactions of the LPWs with bathymetry by O’Kane et al. (2014), such as those studied in idealized numerical models (Barnier 1988). The precise interplay between nonlinearity, bathymetry, and wave propagation is beyond the scope of this paper, but merits further investigation.

The varying propagation speed along the southern pathway is in contrast to the propagation observed along the central (Fig. 4b, topi) and northern (Fig. 4c, top) pathways, where no clear acceleration of the anomalies can be detected. Anomalies on the central pathway take 2–3 years to propagate from the West Australian coast (~100°E) to 80°E. Anomalies on the northern pathway propagate faster still, traversing the basin (from ~110° to ~60°E) in about 1 year. However, a continuum of time scales exists on the northern pathway, with anomalies persisting for as long as 5 years in the western basin, confirming the in-band variance analysis presented in section 1.

To complement the Hovmöller diagrams, we have computed the space–time autospectrum of *ψ*_{g} using Welch’s method at the surface for each of the three pathways (Figs. 4a,b, bottom panels). It is important to note that the spatial coordinate in the autospectrum is the along-path distance (where zero is set to be the western terminus of the path and increases eastward) and that the paths are nonzonal. This is particularly evident in the case of the southern pathway: anomalies propagate northward by ~15° while propagating to the west by ~50° (a ratio of about 0.3). As such, the spatial wavenumber computed from the autospectra is not the zonal wavenumber *k*, but the along-path wavenumber, which we will denote $k^$. We also plot the theoretical dispersion curves from Eq. (1) as thin black lines for the first three baroclinic modes (*m* = 1, 2, 3). In computing the theoretical dispersion curves, we have assumed that $k^\u2248k$ and that *l* ≈ 0, although in the long-wave regime considered here, Eq. (3) suggests that this assumption has a negligible influence on the resulting dispersion curves.

The autospectra show that the dominant frequencies and wavenumbers differ for each propagation pathway. The autospectrum for the southern pathway (Fig. 4a, bottom) has a power distributed among a broad range of frequencies with corresponding periods from about 5 years to about 8 years and across a broad range of wavenumbers (from about 800 to around 2500 km), with most power falling in the left quadrant, indicative of westward propagation. The autospectrum of the central pathway (Fig. 4b, bottom) is similarly broadband, with power distributed across periods from 1 to 6 years, and with wavelengths from about 300 km through to about 1500 km. In contrast to the southern and central pathways, autospectrum of the northern pathway Fig. 4c, bottom) shows power distributed among a range of temporal frequencies (from an equivalent period of about 1 year to around 10 years) but a relatively narrow range of wavenumbers, with a peak wavelength near 3500 km.

In all cases, the autospectra show the greatest degree of similarity to the theoretical dispersion curves corresponding to the first and second baroclinic modes, although the broadband nature of the autospectra in both wavenumber and frequency mean that there is no strong alignment with any particular theoretical dispersion curve. Despite this, we note that the theoretical dispersion curves do provide a reasonable qualitative description of the empirical autospectrum. For example, both the theoretical and empirical autospectra show nondispersive behavior at the wavelengths of interest. Additionally, the quantitative similarity between the theoretical second baroclinic mode estimate of frequency at low wavenumbers provides additional support to conclusions of the feature tracking and CEOF analysis in section 4b that these features are likely a result from a mixture of baroclinic mode 1 and 2 LPWs.

To summarize the results of this section, we find that the anomalous features propagate along distinct preferred pathways. Hovmöller diagrams and spectral analysis of these pathways reveal variability with a range of spatiotemporal scales, with peak periods of around 5 years and peak wavenumbers at long wavelengths (in excess of 1000 km and certainly larger than the local deformation radius). These features show some characteristics of baroclinic LPWs, dominated by modes 1 and 2. Having characterized the low-frequency variability in the South Indian Basin, we now seek to understand its influence on the climatic state.

## 6. How do LPWs influence the state of the climate system

### a. Influence of the LPWs on the oceanic temperature field

The identified mode of variability is primarily a subsurface phenomenon [see Figs. 1c,d(iii)] and the lingering question of whether or not the LPWs are able to influence the climate state through, for example, modification of SSTs, remains to be answered. In this section, we quantify the impact of the LPWs on the near-surface temperature as they propagate along their preferred pathways, and assess their potential impact on the atmosphere.

In Fig. 5, we show a 40-yr-long depth–time slice at two geographic locations on the southern propagation pathways in one point to the east, close to the genesis regions (left panels), and one farther to the west, remote from the genesis region (right panels). Considering the easternmost geographic location on the southern pathway, there is a clear 4–6-yr periodicity in the temperature below about 100-m depth (Fig. 5c), which manifests most notably as excursions in the depth of the 17°C isotherm by as much as 200 m. The deviations in the deep temperature field are strongly correlated with the passage of the LPWs, clearly evident in the anomalous geostrophic velocity [Fig. 5c(i)], with the temperature at a particular depth decreasing (increasing) with anomalous equatorward (poleward) flow, which is consistent with flow in thermal wind balance. However, despite the clear signal of the LPWs at depth, their imprint on the SST is not nearly as clear. Indeed, while there is evidence of 4–6-yr periodicity in the SST at this location, there is an additional high-frequency variability superimposed [Fig. 5a(i)]. Mixed-layer depths (MLDs), computed from the online diagnostics within the ocean model (solid line in Figs. 5b and 5c) are noisier still. The MLDs show negligible correlation with the both the SST and the slow displacements of deep isotherms and appear to have their own, shorter, intrinsic time scale.

Farther to the west along the southern pathway [Fig. 5(ii)], we see behavior similar to that observed to the east: variations of the deep temperature field with a 4–6-yr period and associated anomalous meridional velocities that have some imprint on the SST. However, both the thermocline depth variations and the anomalous velocities have decreased in magnitude and are more episodic (e.g., no anomalies can be detected between years 31 and 36), which suggests that the amplitude of these anomalies decreases as they traverse the basin. Despite this reduction in amplitude, the variance of the SST in the 2–10-yr time bands is approximately the same at both locations, indicating that although the LPWs themselves have lost some coherence and structure, they are still able to influence the ocean surface.

A similar analysis is performed for the central and northern pathways (not shown for brevity). These results mirror those for the southern pathway discussed above, with perturbations of subsurface isopycnals on slow time scales (5–6 years for the central pathway and 1–2 years for the northern pathway) that are concurrent with anomalous geostrophic velocities associated with the propagation of LPWs, and show correlation with the slow modulation of the SST. Like the southern pathway, the central pathway shows decreased variance of subsurface temperatures with distance along the pathway. In contrast, the northern pathway amplitudes intensify in the west of the basin. However, the deep temperature perturbations are also weaker on the northern pathway than either the southern or central pathways, with the variance of the depth of a representative isotherm being less than 50 m, compared to over 150 m for both the southern and central pathways.

To assess time scales of variability present on the LPW pathways, we show the power spectral densities (PSDs) of both the anomalous SST and deep thermocline depth anomalies in Fig. 6 on each pathway. The PSDs are computed at two geographic locations on each LPW pathway using Welch’s method: a point to the east near the Australian continent, and a point to the west near the termination of the pathway. At all locations considered here, both the SST and thermocline depth PSDs show a spectral shape with a relatively flat spectra at low frequencies before transition to a spectral roll off with an approximate slope of *ω*^{−2} at frequencies greater than ~1/2 year (the exact transition frequency being location dependent). This spectral form is strongly reminiscent of an AR(1) process, as discussed in section 2b. However, broad low-frequency peaks are present in both SST and thermocline depth PSDs as well. To assess if these peaks are significant, we test the observed SST PSDs against a null hypothesis of a stochastic AR(1) process (for implementation details, see appendix C). The spectral form of an AR(1) process is identical to that for an ocean mixed layer driven by random, uncorrelated atmospheric forcing (i.e., white noise) with a linear damping, discussed in section 2b [see Eq. (15)]. In essence, we ask if the SST PSDs differ significantly at low frequencies from this simple stochastic model. Note that only the significance of the SST is tested. The results show several of the low-frequency peaks in SST being significant at the 95% level. In particular, the broadband peaks in the SST power spectrum, with accompanying peaks in the thermocline depth PSDs are present at both selected locations on the southern pathway (Figs. 6e,f) and on the western locations on the central and northern pathways (Figs. 6b,d).

The similarity of the deep thermocline and SST PDSs is further investigated by calculating the magnitude squared spectral coherence, the frequency domain equivalent of the squared correlation coefficient (von Storch and Zwiers 1999, 282–288) and the cross-spectral phase, between these two fields (Fig. 7), using Welch’s method, together with the 5%, 50%, and 95% confidence intervals calculated by Monte Carlo methods (see appendix C). In several of the selected locations, we find peaks in the coherence at frequencies between 1 decade^{−1} and 1 yr^{−1}, with the coherence values generally lying between 0.4 and 0.8 in those locations where the peaks exceed the 95% confidence interval. Of the selected locations, only one (the eastern location on the northern pathway; Fig. 7a) has coherence that fails to reach a statistically significant level for time scales between 1 and 10 years, although in several locations (notably the western points on the southern and central pathways, shown in Figs. 7d,f) the coherence is weak. Frequency bands with high coherence generally have either a near-constant or linear increase of phase with increasing frequency. The phase generally falls between 45° and positive 90°, which is indicative of (positive) thermocline depth variations (i.e., deepening) *leading* the SST by a fixed lag. Careful investigation of the depth–time sections, such as that shown in Fig. 5 supports this interpretation: in each 4–6-yr cycle of thermocline depth deepening (shoaling), there is corresponding SST warming (cooling) with the two variables being out of phase by about a quarter of a period (~1–2 years).

To summarize, we have shown that the LPWs that propagate westward across the Indian basin appear as deep-reaching (to around 400-m depth) perturbations in the temperature field. The period and amplitude of these perturbations is dependent on the pathway of the disturbance—being generally stronger and longer period on pathways farther to the south—and the location on each of the pathways. Unsurprisingly, these locations correspond to those locations with the most significant peaks in the power spectra. The disturbances generally remain coherent but reduce in amplitude as they cross the basin. Even though the LPWs manifest most strongly at depth, analysis of the PSDs and spectral coherence of the SST and the deep temperature perturbations indicates some influence of these features on the SST field, which is felt most strongly on the southern pathway, near their genesis regions on western Australian coastline.

### b. Influence of the LPWs on the near-surface ocean

We now investigate the relative influence of the LPWs and atmospheric forcing on SST, including assessment of the potential impact of LPWs on the surface heat fluxes. To accomplish this, we diagnose each term in the MLD temperature budget given by Eq. (7). To compute the geostrophic current velocity *υ*_{g}, we use the horizontal gradient of the 500-db referenced dynamic height anomaly, given by Eq. (18). The MLD, *h*, is taken from the online ocean model diagnostics and is computed as the first level for which the buoyancy change from the surface exceeds 3 × 10^{−4} m s^{−2}. Ekman and surface heat fluxes are taken directly from the atmosphere–ocean model coupler.

Figure 8 shows the in-band variance of the wind forcing [Eq. (11); Figs. 8a,d,g], oceanic forcing [Eq. (12), Figs. 8b,e,h], and the sum of shortwave, longwave, sensible, and latent surface heat flux anomalies [Eq. (10), Figs. 8c,f,i]. The in-band variance of SST (solid contour lines) is also shown in each panel of Fig. 8 as solid lines. Three temporal bands are plotted in Fig. 8: (top) 1–2, (middle) 2–5, and (bottom) 5–10 years.

As is to be expected, the wind forcing (Figs. 8a,d,g) is relatively broadband: similar spatial patterns in variance are observed across all frequency bands. In general, the wind forcing is latitudinally dependent and zonally symmetric, being strong in the tropical region north of 10°S, weak in the subtropical region between 10° and 30°S, and stronger south of about 30°S. The basic zonal symmetry is broken, however, in the 2–5- and 5–10-yr bands, with a region of high variance extending from the central basin to the Australian coast near 20°S. The oceanic forcing (Figs. 8b,e,h) has a similar spatial pattern as the wind forcing but is generally 1 to 2 orders of magnitude stronger. However, in the subtropics, the oceanic forcing organizes into the familiar “beams” associated with the southern and central LPWs pathways, which is particularly clear in the 5–10-yr time band. Areas of high SST variance coincide with areas of high oceanic forcing at low frequency in the subtropics, particularly near 25°S and between 80° and 100°E. The in-band variance of the surface heat flux (Figs. 8c,f,i) shows higher variance in tropical regions at frequencies less than 5 years. In the subtropics, the high variance regions align well with regions of high SST variance, particularly in the 2–5- and 5–10-yr time bands. The variability seen in the 5–10-yr band coincides particularly well with the variability in the SST and shows remarkable similarity with the patterns of oceanic forcing south of about 20°S. These two fields are most similar in regions corresponding to the southern and central LPW pathways.

The similar spatial patterns in the oceanic forcing and surface heat fluxes in 5–10-yr time bands along certain LPW pathways strongly hints that these two fields are both related to the passage of these propagating features. To investigate further, we plot 40-yr Hovmöller diagrams for the wind forcing, oceanic forcing, and surface heat flux anomalies along the southern LPW pathway in Fig. 9, which shows the greatest apparent similarity between the surface heat flux and oceanic forcing. We have low-pass filtered each quantity using a finite-impulse-response Blackman filter with a cutoff period of 1 year (although the oceanic forcing in Fig. 9b is essentially unmodified by the filtering operation). The wind and oceanic forcing are converted into equivalent units of watts per square meter by multiplication by *ρ*_{0}*C*_{p}*H*, where *H* is the time mean mixed-layer depth along each pathway. For all variables, westward propagation can be seen with a period of between 4 and 6 years, consistent with the LPWs discussed in section 5, although we note that the wind forcing anomalies are both weaker than the oceanic forcing anomalies and appear to propagate more rapidly to the west than either the oceanic forcing anomalies or the surface heat flux anomalies. Surface heat flux anomalies, in contrast, propagate westward with an almost identical phase speed to the oceanic forcing anomalies, and have similar magnitudes, although we note that close study of Fig. 9 reveals that the fields are out of phase by approximately 1–2 years. Heat fluxes are dominated by latent heat fluxes, which we find to be approximately twice as strong as either sensible heat fluxes or longwave radiation fluxes, although the individual components of the fluxes most often act in concert rather than opposition. Anomalous geostrophic heat flux divergence and entrainment terms contribute equally to the oceanic forcing (although entrainment is somewhat more episodic). Other terms in Eq. (12) are negligible.

Heat flux anomalies associated with the oceanic forcing have a peak magnitude of ~30 W m^{−2}, while anomalous fluxes with magnitudes greater than 20 W m^{−2} may persist for as long as 2 years in certain locations. Additionally, we note that there are periods where the heat flux anomaly appears to be controlled by the oceanic forcing. For example, between model years 25 and 35, the heat-flux signature strongly resembles that of the oceanic forcing, albeit with a significant phase shift. In contrast, the wind forcing of the SST in this period does not show a strong signal.

The similarity of surface heat flux and oceanic forcing Hovmöller diagrams shown in Fig. 9 provides further evidence of a substantial influence of slow subsurface oceanic variability on the surface ocean, at least along the southern pathway. To quantify the strength of this relationship on all the LPW pathways considered here, we plot in Fig. 10 the magnitude squared spectral coherence and cross-spectral phase between the total surface heat flux and the deep thermocline anomalies, at two locations (east and west) on each LPW pathway, as described in section 5. We also include estimates of the 5%, 50%, and 95% confidence intervals (see appendix C). Statistically significant correlations greater than 0.5 are found at frequencies between 5 and 9 years at the eastern locations on both the central and southern pathways, as well as at the western location on the northern pathway. In each of these locations, the phase remains within a relatively narrow band of values, being generally 90° to 180° (i.e., one-quarter to one-half of a cycle) out of phase, with (negative) thermocline depth anomalies leading the surface heat fluxes. In all other locations, and at higher frequencies, the coherence either fails to be statistically significant or is small (typically less than 0.2).

The regions with high coherence at low frequencies are unsurprisingly also those with relatively high power at low frequencies in the SST PSDs (Fig. 6) and strong coherence between the SST and the deep thermocline anomalies (Fig. 7). The coherence analysis undertaken in this section and in section 6, together with the in-band variance presented in Fig. 8, is strongly suggestive of a relationship between the deep thermocline anomalies associated with the passage of LPWs and anomalous surface heat fluxes. At locations remote from the LPW genesis, surface heat flux anomalies in excess of 20 W m^{−2} for periods of 1–2 years.

## 7. Discussion

### a. Can LPWs induce a response in the atmospheric circulation?

In a coupled climate model we have demonstrated that simulated LPWs traverse the Indian Ocean Basin and are capable of influencing the near-surface ocean. Indeed, we have shown that persistent heat fluxes between the ocean and atmosphere on the order of 20 W m^{−2} are induced by the passage of LPWs along the southern pathway. For comparison, the total climatological heat flux in our model in the subtropical south Indian Ocean ranges from ~+50 W m^{−2} in summer to ~−150 W m^{−2} in winter. As such, the anomalous heat flux induced by the LPWs is around 50% of the summer climatological values, and about 20% of the winter climatological values. These values are significantly smaller than those used in idealized studies of the atmosphere driven by imposed SSTs, where the imposed anomalous fluxes are typically 500–1000 W m^{−2} (e.g., Kaspi and Schneider 2011). Are the relatively small heat fluxes associated with the Indian Ocean LPWs large enough to generate an organized atmospheric response? Here, we provide evidence that suggests LPWs are indeed capable of forcing an atmospheric response by exploiting our fully coupled atmosphere–ocean model.

We plot the 40-yr time series of the anomalous atmospheric vertical (pressure) velocity *ω* = *Dp*/*Dt* in Fig. 11a at the eastern location on the southern LPW pathway, low-passed filtered with a cutoff period of 1 year. We note that low-pass filtering reduces the atmospheric variance by about 80%, as the region is dominated by high-frequency synoptic features. Clear low-frequency variability is present in both the wind divergence and *ω* fields, each with a period of between 2 and 6 years.

To link the near-surface atmospheric response to the surface ocean, we show in Fig. 11b, the anomalous total surface heat flux. The near-surface atmospheric variability is strongly correlated with variability in the surface heat flux, in which periods of anomalous positive heat flux (the sign convention of our model is that positive heat flux implies a warming of the ocean) tend to coincide with positive values of *ω*, indicative of downward motion. To relate the atmospheric motion to the ocean dynamics, we plot in Fig. 11c the oceanic conservative temperature from the surface to 400-m depth. As discussed in section 6, there is a clear correlation between the between the depth of the thermocline (shown as the dashed line in Fig. 11c) and the surface heat flux into the ocean: as the thermocline shallows, colder water is brought to the ocean surface and the anomalous surface heat flux becomes more positive (indicative of the ocean gaining heat from the atmosphere). As the thermocline deepens, the surface ocean warms and, in turn, the heat flux into the ocean becomes more negative (indicative of the ocean losing heat to the atmosphere). The surface heat flux lags the thermocline depth anomalies by around 1 year.

To further investigate the correlation between oceanic and atmospheric variability, we compute the magnitude squared coherence between the deep thermocline anomalies and the atmospheric vertical velocity *ω* (Fig. 12) at the eastern and western locations on the LPWs pathways discussed previously, at all atmospheric levels between the surface and 250 hPa. Stippling indicates correlations significant at the 95% level (determined following the procedure outlined in appendix C, applied to each level separately). Significant correlations between the atmospheric vertical velocity and the deep thermocline anomalies greater than 0.5 are found in frequency bands corresponding to periods between 5 and 10 years, which falls within the time bands expected for LPWs. Significant correlations extend to a range of atmospheric heights, to 800 hPa in the case of the eastern location on the southern pathway, to 500 hPa on the central pathway, to 300 hPa on the northern pathway. In contrast to eastern locations, the western locations do not show any appreciable correlations, even in regions where there is strong coupling between the subsurface and surface ocean (such as on the northern pathway). At present, we can only speculate on the reasons for this diversity in behavior. However, we note that Fig. 8 shows that the oceanic forcing of the SST is weaker west of 80°E at time scales greater than 2 years. As such, although the deep thermocline anomalies may be well correlated with surface heat fluxes in some of these regions, the response of the surface ocean may simply be too weak to generate an organized atmospheric response.

From the brief analysis presented here, we have shown that—at a few geographic locations—low-frequency variability of the deep thermocline is correlated with low-frequency variability in the atmosphere. The time scales of the correlation correspond well to those associated with the propagating LPWs. From this analysis, we tentatively conclude that the low-frequency variability associated with the LPWs may be capable of influencing the atmosphere in locations where oceanic influence on the surface heat flux is strong. Understanding the detailed dynamics of the ocean–atmosphere interaction is beyond the scope of this work, but will be explored in future work.

### b. Potential atmosphere–ocean feedbacks

Using observational and reanalysis products, White (2000) found significant covariance between the surface wind field and westward-propagating SST anomalies in the Indian Ocean. He argued that the observed covariance was evidence of a significant ocean–atmosphere feedback, that partially explained the fact that the observed SST anomalies propagated significantly slower than the theoretical first-baroclinic-mode phase speed would suggest. Having shown that the slow ocean circulation is capable of influencing the local atmospheric circulation at several locations on the eastern side of the basin, it is natural to ask if ocean–atmosphere feedbacks might play a role in the maintenance and propagation of the modeled LPWs.

Following White (2000), we plot Hovmöller diagrams of the meridional 10-m wind anomalies in Figs. 13a–c along each of the three LPW pathways for a 40-yr time period. The wind anomalies have been low-passed filtered, with a cutoff period of 1 year. While there is evidence of low-frequency variability and coherent propagation of the 10-m wind anomalies, both eastward and westward propagation are found. Additionally, the propagation of the surface wind anomalies is substantially faster than that observed in the ocean, regardless of the propagation pathway and the direction of travel. For example, the basin crossing time is generally less than one year on the southern pathway, compared to almost 10 years for the oceanic propagating features.

To complement the Hovmöller diagrams, we compute the autospectrum of the 10-m meridional wind anomalies, shown in Figs. 13d–f. The 10-m meridional wind autospectra show little similarity to those of the ocean dynamic height, being broadband in frequency, but concentrated at low wavenumbers. The autospectra also show substantial power in the positive wavenumbers, which corresponds to eastward propagation.

Analysis of the propagation of the 10-m meridional wind anomalies along the LPW pathways and their associated autospectra reveals a lack of clear covariance with the oceanic LPWs and strong contrast in their propagation speed and direction. As such, we find no clear evidence of the ocean–atmosphere feedback studied by White (2000) in this model. We hypothesize that while the low-frequency variability in the ocean may force the atmosphere in some geographic locations, the excited atmospheric circulation does not necessarily feed back on the oceanic circulation.

### c. Implications for predictability

The LPWs discussed in this study are long-lived features that remain coherent despite propagating long distances. As we have shown, the LPWs are also capable of substantially modifying the local environmental conditions. As such, if these structures were predictable in some sense, they could be used as a route to improve predictability on seasonal or even decadal time scales. However, as noted by Wolfe et al. (2017), despite the LPWs having a fairly regular period and phase speed, they have highly irregular amplitudes and phases. Wolfe et al. (2017) investigated the predictability of this system by running a parallel model experiment: two forced ocean models that differed only by a small (order 10^{−7}) perturbation applied to their initial conditions were integrated forward in time for 60 years and the resulting divergence of the trajectories of the two experiments projected onto the mode of Indian Ocean variability discussed in this paper. They found that the model trajectories diverged at an exponential rate, with an *e*-folding time scale around 2 years. Indeed, the authors noted that the intrinsic variability associated with the propagating LPWs in the Indian Ocean tends to overwhelm the forced response on time scales longer than a few years.

While Wolfe et al. (2017) showed exponential divergence in their perturbation experiments, we find their conclusions somewhat too pessimistic. In particular, we note that their perturbation experiment is likely to excite small length scale and short time scale instabilities associated with the first finite size Lyapunov exponent of the system. However, the predictability in geophysical systems typically lies in larger-scale, slower-time-scale phenomena (Boffetta et al. 1998). More modern ensemble generation techniques, such as bred vectors (Toth and Kalnay 1997), allow one to suppress the growth of small-length-scale, short-time instabilities relative to the large-space-scale, slow-time-scale dynamics. Such an approach was used by Hoffman et al. (2009) to study the dynamics of equatorial waves and by O’Kane et al. (2011) to study the potential predictability of the East Australian Current.

The predictability of the Indian Ocean LPWs is beyond the scope of this study. However, we note that as relatively large-scale structures that remain coherent over long periods of time, they should provide an ideal target for medium-term climate forecasting studies. Full predictability studies in this phenomenon, including the application of bred vector techniques, is the subject of ongoing work.

## 8. Conclusions

In this study, a fully coupled climate model has been used to study low-frequency variability in the Indian Ocean. On time scales of 2–10 years, the variability is dominated by propagating features that follow one of several distinct pathways that cross the Indian Ocean from east to west. The features that propagate along these pathways remain relatively coherent and have dispersion characteristics consistent with long planetary waves (LPWs). Although the modeled waves are broadband in both frequency and wavenumber, with substantial modulation of amplitude and phase, analysis of the propagation speed of the tracked features, supplemented by CEOF analysis, leads us to conclude that these features are likely a mixture of baroclinic mode 1 and 2 LPWs, which is in agreement with recent work by O’Kane et al. (2014) and Wolfe et al. (2017).

The LPWs manifest primary as large excursions in the depth of the thermocline by as much as 200 m of a representative isotherm. Although the variability manifests primarily at depth, in certain geographic locations the LPWs can significantly impact the SST and the surface buoyancy fluxes both into and out of the atmosphere. By using a one-dimensional mixed-layer temperature budget, we are able to separate the contributions of the fast wind forcing and the slow ocean forcing. We find that, on time scales longer than about 2 years, oceanic forcing dominates variance in the SSTs over much of the Indian Ocean, particularly on the eastern side of the basin. The deep thermocline variations are strongly correlated in certain locations with anomalous surface heat fluxes, with amplitudes of ±20 W m^{−2} that are capable of being sustained for periods of as long as 2 years. We have presented evidence that this anomalous surface flux is, in turn, capable of generating a coherent atmospheric response, at least in the lower atmosphere, although there does not appear to be a strong feedback to the ocean. The mode of near-term variability and its influence on the surface ocean and near-surface atmosphere is summarized schematically in Fig. 14.

The mode of variability studied here resembles the variability described by O’Kane et al. (2016) and Wolfe et al. (2017), who both considered an ocean only model driven by repeat-year atmospheric forcing: a single “normal” year (i.e., ENSO neutral) that repeated year after year. This forcing was selected specifically to isolate the intrinsic variability in the ocean, as models forced with repeat-year forcing will not generate ENSO-like signals that are known to have a significant impact on the Indian Ocean. The fact that both O’Kane et al. (2016) and Wolfe et al. (2017) were able to generate a similar mode of variability using models with no interannual variability in forcing provides evidence that the features studied in this paper results from an intrinsic ocean variability, and are unrelated to the projection of variability arising in the Pacific influencing the Indian Ocean Basin. However, as ENSO variability projects strongly into the Indian basin in similar time bands (4–5 years) to that associated with the LPWs, it is likely that ENSO could modulate, suppress, or otherwise interact with the generation of LPWs on the eastern boundary of the basin. This effect is likely to be strongest in the tropical band (i.e., along the northern pathway), but could also be communicated to the central and southern pathways by the excitation of Kelvin waves along the western coast of Australia.

In this study, we have not investigated in detail the dynamics of the generation of LPWs or the dynamics of their propagation. Wolfe et al. (2017) present energetics and linear-stability arguments that suggest that they are the result of quasigeostrophic boundary instabilities. However, other plausible mechanisms, such as eastern boundary upwelling and excitation by anomalous variations in the Leeuwin Current and the Indonesian Throughflow, should also be investigated. Additionally, the roles of topography, background flow, and nonlinearity on LPWs have not been addressed in this study, yet they could play crucial roles in governing their propagation characteristics and their capacity to support long-range ocean teleconnections.

We have also hinted at a coherent atmospheric response to the LPWs. White (2000) and White et al. (2004) noted using satellite observations and atmospheric reanalysis that LPWs can couple with the atmosphere, driving a response in the near-surface winds that, in turn, can modify the LPWs themselves. However, the complete coupled response, including the large-scale atmospheric response to the localized oceanic forcing and any feedbacks between the two systems, has also not been addressed in this study.

Finally, the analysis here has been conducted in a relatively low-resolution climate model that is incapable of representing important features of the climate systems, such as mesoscale eddies. Robust detection and characterization in either observations or high-resolution numerical models could shed further light on important mechanisms governing the near-term variability in the ocean and their potential predictability. These questions will be addressed in future work.

## Acknowledgments

The authors were supported by the Australian Commonwealth Scientific and Industrial Research Organisation Decadal Forecasting Project (research.csiro.au/dfp). The authors would also like to acknowledge NCI for providing the computational resources. The comments of three anonymous reviewers greatly improved the manuscript.

### APPENDIX A

#### Feature Detection and Tracking

The feature detection and tracking used to define the propagation pathways of the LPWs is based on the mesoscale eddy detection and tracking algorithm of Chelton et al. (2011). The algorithm proceeds by searching for simply connected regions enclosed by a single contour of anomalous dynamic height. The advantage of this approach, compared to more common approaches such as the use of two-dimensional Radon transforms (e.g., Maharaj et al. 2009), is that our approach allows the detected features to propagate in two dimensions. Importantly, this approach does not require a prescribed dynamic height threshold contour to be defined: to identify a feature, all that is required is a region enclosed by single dynamic height anomaly contour, subject to a few different criteria used to filter features that are too large, too small, or too weak.

The major steps of the feature detection algorithm are as follows:

Starting with a snapshot of anomalous dynamic height at the ocean surface,

*ψ*_{g}(*x*,*y*,*p*= 0,*t*), we identify all simply connected regions that exceed (in the case of anticyclones) or are smaller than (in the case of cyclones) a given critical dynamic height threshold, such that $\psi g\u2032>\psi gcrit$ (for anticyclones) or $\psi g\u2032<\psi gcrit$ (for cyclones) (see Fig. A1a).Each identified region is evaluated against a set of criteria (listed below) to ensure that it is neither too large nor too small, and to check if it contains at least one local extremum (Fig. A1b). Those regions that do not meet the specified criteria are rejected.

The center of mass (or centroid) of each of the remaining regions determines the latitude and longitude of the identified feature (Fig. A1c). All pixels contained within the $\psi gcrit$ contour are masked.

Steps 1 to 3 above are repeated for a new, smaller absolute value of $\psi gcrit$. To find the largest, outermost contour enclosing the feature we test values of $\psi gcrit$ from 100 to −100 cm in increments of 5 cm. Anticyclones and cyclones are detected separately.

As mentioned above, regions identified as being enclosed by a dynamic height anomaly are filtered based on the following criteria:

has a magnitude that exceeds 1 cm;

contains at least one local maximum;

covers an area greater than 4 grid boxes, but no more than 60 grid boxes.

The first two criteria are used to separate weak and generally noisy regions from coherent LPWs. The third criterion removes regions that have length scales too small or too large to be LPWs (although in practice, regions are rarely excluded for being too small). The exact values selected for these parameters were informed by the theory of linear planetary waves, but have been obtained primarily through trial and error.

Tracks of coherent features are formed by comparing the locations of detected features between successive monthly time steps: two features at time step *n* and *n* + 1 of the same sign (either cyclonic or anticyclonic) are declared to be part of the same “feature” if they are found within an ellipse with an east–west axis length of 600 km, and a north–south axis length of 500 km. As above, the parameter values selected were informed by theory, but were selected primarily by trial and error.

The values of the hyperparameters in the algorithm above are modified from the values used in studies of mesoscale eddies. To determine the robustness of the algorithm, we have run a series of sensitivity tests where the values of the threshold amplitude, minimum peak amplitude, area thresholds, and the axis lengths of the tracking ellipse were all varied systematically. We find that the broad-scale results are largely insensitive to the values of the hyperparameters provided that each lies within a rather large physically plausible envelope. The exceptions to this statement are the value of the threshold amplitude and the tracking ellipse major (east–west) axis. We have found that these parameters must be selected carefully; otherwise, the methodology will either fail to identify any LPWs in the domain or will not track the features over successive time steps.

The implementation was taken from the open source software described in Oliver et al. (2015), which can be obtained from the author’s Github repository (https://github.com/ecjoliver/eddyTracking).

### APPENDIX B

#### Derivation of the Anomalous Mixed-Layer Temperature Budget

Here, we derive the complete anomalous mixed-layer temperature budget [Eq. (7)]. Although similar budgets have been utilized by previous authors, many derivations are somewhat ad hoc, and neglect certain terms in the expression without complete justification.

We begin with the mixed-layer temperature budget [Eq. (5)], rewritten here:

Each dependent variable is expanded as a seasonally varying time mean $(\u22c5)\xaf$, and a perturbation (·)′ (e.g., $\Theta =\Theta \xaf+\Theta \u2032$, $h=h\xaf+h\u2032$, …). We note the identities:

After expansion, making use of identities (B1) and (B2) and the incompressible continuity equation ∇_{xy} · **u** + (∂*w*/∂*z*) = 0, we take the time mean of Eq. (5) to obtain an expression for the seasonally varying mixed-layer temperature:

To derive an expression for the anomalous mixed-layer temperature, we subtract Eq. (B3) from the expanded Eq. (5), which gives, after some algebra,

Terms describing covariances between two perturbation variables with the form $a\u2032b\u2032\xaf$ have maximum periods of 1 year due to the averaging operation. As such, we ignore them as they are irrelevant for the low-frequency variability studied in this paper. Ignoring these terms and dividing by $h\xaf$ gives Eq. (7).

### APPENDIX C

#### Confidence Intervals for the Power Spectral Densities and Spectral Coherence

In this appendix, we describe the method for calculating the confidence intervals used to determine statistical significance for both power spectral densities (PSDs) and the magnitude squared coherence. In both cases, physically intuitive tests using Monte Carlo methods are favored over formal hypothesis tests, due to their relatively straightforward interpretation.

Beginning with the PSDs, we remind the reader that the goal of this paper has been to determine the influence of interannual planetary waves on the near-surface atmosphere and ocean. If these phenomena are to have an influence, they should manifest at low frequencies significantly different to an ocean or atmosphere driven by high-frequency forcing. Following Frankignoul (1985), we assume as a null hypothesis a first-order autoregressive process, which has the PSD

where *λ* is the damping parameter and *P*_{0} is the PSD of a white stochastic process (constant at all frequencies) of (as yet) undetermined magnitude. This PSD is identical to the spectrum of a simple mixed layer driven by atmospheric forcing that has a white noise spectrum (Frankignoul 1985).

For a given PSD, the parameters *λ* and *ω* are estimated in the frequency domain by using a de-biased Whittle maximum likelihood (Sykulski et al. 2019). We then generate a random sequences with a PSDs that follow Eq. (C1) by first generating a white stochastic sequence with unit variance. We then Fourier transform this sequence into the frequency domain using the fast Fourier transform, which results in a PSD with an average power (across all frequencies) of 1, but random phase angles. We multiply both the real and imaginary Fourier coefficients by the square root of the fitted PSD (which corresponds to the signal amplitude) and then apply the inverse fast Fourier transform to generate a random AR(1) time series with the desired autocorrelation structure

This process is repeated 10 000 times, and the PSD of each of the realizations is calculated using Welch’s method. From these random PSDs, we are able to calculate the confidence intervals by determining the 5th, 50th, and 95th percentiles.

Confidence intervals are computed for the magnitude square coherence spectra in a similar manner. First, if we consider the coherence of two time series, *x*(*t*) and *y*(*t*), we test their calculated coherence against the null hypothesis of two uncorrelated time series with the same autocorrelation structure as *x* and *y*. To do this, we generate 10 000 random sequences as before. However, instead of using the PSD fitted to Eq. (C1), we use the true PSDs, smoothed in the frequency domain using a boxcar filter spanning three frequency bins. We then compute the magnitude squared coherence between the random sequence and the original variables. This procedure is repeated 10 000 times to determine the 5%, 50%, and 95% confidence intervals as before.

## REFERENCES

*J. Phys. Oceanogr.*

*Aust. Meteor. Oceanogr. J.*

*Phys. Chem. Earth*

*Physica D*

*Sociol. Methods Res.*

*Science*

*Prog. Oceanogr.*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*J. Climate*

*J. Geophys. Res.*

*Rev. Geophys.*

*Tellus*

*J. Phys. Oceanogr.*

*Wiley Interdiscip. Rev.: Climate Change*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*Atmosphere–Ocean Dynamics*. Academic Press, 662 pp.

*J. Climate*

*Climate Dyn.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Bull. Amer. Meteor. Soc.*

*Tellus*

*Geophys. Res. Lett.*

*J. Atmos. Sci.*

*J. Geophys. Res.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Ocean Dyn.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Ocean Modell.*

*Dyn. Atmos. Oceans*

*J. Geophys. Res. Oceans*

*Geophys. Res. Lett.*

*J. Geophys. Res.*

*Control Systems Engineering*. 6th ed. John Wiley & Sons, 926 pp.

*Ocean Modell.*

*J. Geophys. Res. Oceans*

*J. Geophys. Res. Oceans*

*J. Climate*

*J. Geophys. Res. Oceans*

*Numerical Recipes 3rd Edition: The Art of Scientific Computing*. 3rd ed. Cambridge University Press, 1235 pp.

*J. Phys. Oceanogr.*

*Rev. Geophys.*

*J. Geophys. Res.*

*Geophys. Res. Lett.*

*J. Mar. Res.*

*J. Phys. Oceanogr.*

*Biometrika*

*Mon. Wea. Rev.*

*J. Geophys. Res. Oceans*

*Statistical Analysis in Climate Research*. Cambridge University Press, 484 pp., https://doi.org/10.1017/CBO9780511612336.

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Climate*