## Abstract

Lagrangian statistics of the surface circulation in the Indian Ocean (IO) are investigated using drifter observations during 1985–2013. The methodology isolates the influence of low-frequency variations and horizontal shear of mean flow. The estimated Lagrangian statistics are spatially inhomogeneous and anisotropic over the IO basin, with values of ~6–85 × 10^{7} cm^{2} s^{−1} for diffusivity, ~2–7 days for integral time scale, and ~33–223 km for length scale. Large diffusivities (>20 × 10^{7} cm^{2} s^{−1}) occur in the central-eastern equatorial IO and the eastern African coast. Small diffusivities (~6–8 × 10^{7} cm^{2} s^{−1}) appear in the subtropical gyre of the southern IO and the southeastern Arabian Sea. The equatorial IO has the largest zonal diffusivity (~85 × 10^{7} cm^{2} s^{−1}), corresponding to the largest time scale (~7 days) and length scale (~223 km), while the eastern coast of Somalia has the largest meridional diffusivity (~31 × 10^{7} cm^{2} s^{−1}). The minor component of the Lagrangian length scale is approximately equal to the first baroclinic Rossby radius (*R*_{1}) at midlatitudes (*R*_{1} ~ 30–50 km), while the major component equals *R*_{1} in the equatorial region (*R*_{1} > 80 km). The periods of the energetic eddy-containing bands in the IO in Lagrangian spectra range from several days to a couple of months, where anticyclones dominate. A significant result is that the drifter-derived diffusivities asymptote to constant values in relatively short time lags (~10 days) for some subregions of the IO if they are correctly calculated. This is an important contribution to the ongoing debate regarding drifter-based diffusivity estimates with relatively short Lagrangian velocity time series versus tracer-based estimates.

## 1. Introduction

Lagrangian statistics of the ocean circulation, including diffusivity, integral time scale, and length scale, determine the dynamic nature of diffusive processes in the oceans. The Lagrangian integral time or length scales represent a “memory” scale following the particles over which the Lagrangian velocities are correlated, that is, over such a time or length scale, a particle “remembers” its past motion state and shows a strong autocorrelation. Long after (well beyond) that time scale (length scale), the particle “forgets” its past state and shows statistically uncorrelated and independent motion (LaCasce 2008). In other words, the diffusion is in a ballistic regime in a time period much shorter than the time scale, while it is in a random-walk regime after times much longer than the time scale (Lumpkin and Flament 2001). The Lagrangian integral time scale and length scale are usually assumed to be related to each other through the characteristic speed, that is, the time scale is the time a particle takes to travel the length scale at the characteristic speed of the turbulent flow (e.g., Lumpkin et al. 2002). These scales are thus essential for determining the values of diffusivity coefficients used in ocean circulation simulations, especially in coarse-resolution models that do not explicitly resolve eddy transport processes (e.g., Bauer et al. 1998; Stammer 1998).

The diffusivity of the ocean circulation can be estimated from Eulerian data (e.g., Bryden and Heath 1985; Ferrari and Nikurashin 2010; Holloway 1986), from tracer data (e.g., Armi and Stommel 1983; Ferrari and Polzin 2005; Jenkins 1998; Klocker et al. 2012b; Zika et al. 2010), and from Lagrangian (particle) data (e.g., de Verdiere 1983; Davis 1991; Krauss and Böning 1987; Lumpkin et al. 2002; Oh et al. 2000; Taylor 1922). However, there is still much recent controversy regarding the discrepancy between different estimates. For instance, Sundermeyer and Price (1998) showed their Lagrangian estimates of diffusivity were 2–6 times larger than the tracer-based estimates using an idealized ocean mode as well as synthetic floats and tracers. In contrast, Riha and Eden (2011) found their Lagrangian and Eulerian diffusivity estimates were comparable, with the Lagrangian ones being slightly smaller. More recently, Klocker et al. (2012b) examined further the relationship between different diffusivity estimates using synthetic particles and tracers advected by a velocity field derived from sea surface height measurements in the South Pacific and demonstrated that particle- and tracer-based estimates of eddy diffusivity are equivalent, assuming sufficient integration times are permitted in the Lagrangian approach. However, given that obtaining sufficient sampling of deployed tracers are costly, and simulated currents from altimetry are generally too low in energy compared to in situ observations, the float-/drifter-based approach is still the most desirable (Klocker et al. 2012b).

Using observations by satellite-tracked surface drifters, many studies have investigated the Lagrangian statistics of near-surface circulation over different ocean basins (e.g., Bauer et al. 2002; Bograd et al. 1999; de Verdiere 1983; Krauss and Böning 1987; Lumpkin and Flament 2001; Lumpkin et al. 2002; Sallée et al. 2008; Veneziani et al. 2004; Zhurbas and Oh 2003, 2004) as well as marginal seas (Centurioni et al. 2004, 2009; Oh et al. 2000; Poulain 2001; Qian et al. 2013; Rinaldi et al. 2010). Many approaches of deriving the Lagrangian statistics of ocean currents from drifter observations are based on the classical framework of Taylor (1922), that is, an integration of the autocorrelation function of drifter-observed velocity. However, such an integration of the autocorrelation function can be quite sensitive to the asymptotic behavior of the autocorrelation at large time lags (Klocker et al. 2012b; Lumpkin et al. 2002). Griffa et al. (1995) tried to avoid this problem by fitting the autocorrelation with a known-shape curve prescribed by some parameters and then only estimating the parameters. To make Taylor’s approach applicable to inhomogeneous mean flows, removal of the mean flows is necessary and should be carefully designed (Davis 1987, 1991; Oh et al. 2000). Bauer et al. (1998) used a bicubic spline interpolation scheme to fit the mean flows and then removed them to minimize the mean shear-induced dispersion. Oh et al. (2000) showed that a principal axis decomposition of the traditional diffusivity tensor could transfer the shear-induced dispersion to the major component, leaving the minor component of diffusivity a good approximation of the asymptotic value. In addition to the spatial variation, the temporal variation of mean flows may have significant impacts on the Lagrangian statistics (Qian et al. 2014; Zhurbas et al. 2014). Lumpkin and Johnson (2013) used the Gauss–Markov (GM) method to estimate the mean, temporal variability, and spatial shear of flows and then removed these components from the total observed velocity to obtain the residual velocity. Qian et al. (2014) further examine the GM method in the diffusivity estimate with synthetic drifters in both the idealized and satellite-based flows in the Indian Ocean (IO) basin. Their results demonstrate that a temporally and spatially continuous fit through the GM estimator is very efficient in isolating the effects of both the shear dispersion and seasonal variation of mean flows, resulting in the best estimate to the “true” diffusivity, especially in regions where strong signals of seasonal (both annual and semiannual) cycles exist, such as the eastern coast of Somalia and the equatorial IO subject to the Indian monsoons. Therefore, this method is particularly suitable for the IO where significant seasonal variability and large horizontal shear of the surface currents occur as indicated in the companion paper (Peng et al. 2015, hereinafter PQL1) and previous studies (e.g., Schott and McCreary 2001; Schott et al. 2009).

Although the Eulerian statistics of the near-surface circulation in the IO have been investigated in some studies using observations from ship reports or surface drifters (e.g., Beal et al. 2013; Cutler and Swallow 1984; Molinari et al. 1990; Shenoi et al. 1999; Zheng et al. 2012), few studies have investigated the Lagrangian statistics of the near-surface circulation in the IO basin. One of the reasons for this could be that the relatively sparse distribution and short span of the drifter observations in the IO was until quite recently insufficient for estimating Lagrangian statistics. Another reason could be that large horizontal shear of near-surface currents and significant seasonal variability occur in the IO that have been shown to influence greatly the estimate of Lagrangian statistics (Bauer et al. 1998; Qian et al. 2014; Zhurbas et al. 2014). Quite recently, Chiswell (2013) provided an estimate of diffusivity for the 1000-m level and the surface circulation over the IO using the half growth rate of dispersion that is based on particle displacements. However, as it is not easy to decompose displacement temporally (displacement is defined between one point and another reference point), Chiswell (2013) did not take into account the seasonal variability, and his estimate over the surface equatorial regions of both the Indian and Pacific Oceans reached ~200 × 10^{7} cm^{2} s^{−1}. This estimate is comparable to that obtained by the traditional binning method without excluding seasonal cycles with synthetic drifters as shown in Qian et al. (2014), which is about 3–4 times larger than the true value for the synthetic drifters (~50 × 10^{7} cm^{2} s^{−1}). Zhurbas et al. (2014) mapped the lateral diffusivity over the globe with emphasis on the IO using Eq. (4) through a principal axis decomposition of the traditional diffusivity tensor (Oh et al. 2000) plus a seasonal binning method (12 temporal bins) to reduce the effects of both inhomogeneous and nonstationary mean flow. However, Zhurbas et al. (2014) did not explicitly remove the shear inside bins, and the monthly density of observations in a bin is far lower than the overall density used in the GM approach. By adopting a principal axis decomposition, the shear dispersion is effectively transferred to the major component, leading to a biased major component of diffusivity. Most recently, although Qian et al. (2014) examine the effectiveness of the GM method in isolating the effects of both the shear dispersion and the seasonal variation of mean flows on the diffusivity estimate in the IO basin, the data they used are synthetic drifters instead of real ones with either idealized or satellite-based background mean flows. Therefore, a reanalysis of the diffusivity tensor as well as its asymptotic behavior with time in the IO by applying the GM method to real drifter data is imperative and will provide new insights from a different perspective than that of Zhurbas et al. (2014). As described in the companion paper PQL1, at present we have accumulated sufficient drifter observations in the IO to resolve annual and semiannual fluctuations and lateral shear in most of the basin. This dense dataset, facilitated with an advanced technique developed in recent years that can effectively remove the spatial shear (inhomogeneity) and temporal variability (nonstationarity) of near-surface currents (Lumpkin and Johnson 2013; Qian et al. 2014), sheds new light on the feasibility of such an investigation and thus motivates the present study.

The rest of this paper is organized as follows: Section 2 gives a description of the drifter data and methods. Section 3 displays the probability density function of the residual velocities. The estimated diffusivities and integral scales of near-surface currents in different regions are presented in section 4. Section 5 discusses the Lagrangian spectra of eddy velocities. A summary is given in the final section.

## 2. Data and methods

The Global Drifter Program (GDP) drifter dataset from 1985 to the most recent data (June 2013), provided by the Atlantic Oceanographic and Meteorological Laboratory/National Oceanic and Atmospheric Administration (AOML/NOAA) (www.aoml.noaa.gov/phod/dac/dacdata.php), is used in this study. This dataset includes quality-controlled positions and velocities that are interpolated to uniform 6-h intervals. A detailed data description, including the data processing as well as data’s spatial/temporal distribution, is given in PQL1.

The drifter data are first binned into 2° latitude/longitude grids, which are chosen as a compromise between sample numbers within a bin and the resolution describing the mean flow. Then the GM decomposition (Lumpkin and Johnson 2013) is applied to each bin to obtain the residual velocity by removing the spatially varying mean and temporally varying seasonal variability of the velocities within a bin. After obtaining , the Lagrangian statistics, including Lagrangian diffusivity matrix as well as the Lagrangian integral time scale and length scale , can be computed as (Davis 1991; Poulain 2001)

using the pseudotrack method (e.g., Swenson and Niiler 1996), where is the Lagrangian ensemble average over time and space at time lag before/after the particles are located in a given bin, and the velocity covariance matrix.

It is worth noting that an alternative way to compute the diffusivity tensor is based on the velocity and displacement residuals:

where is the departure from Lagrangian ensemble-mean displacement (e.g., Oh et al. 2000; Swenson and Niiler 1996; Zhurbas et al. 2014). Although it is theoretically equivalent to Eq. (1), both methods in practice could result in different estimates as the averaging is applied at different stages (LaCasce et al. 2014). We choose Eq. (1) in the present study because it relates diffusivity only to the velocity residual rather than both the velocity and displacement residuals. Seasonal cycles of velocity can then be removed through the GM decomposition, whereas those of the displacement cannot be removed directly. Therefore, effects of seasonal cycles on diffusivity estimation could be isolated more completely using Eq. (1). This urges a comparison of the estimates in the present study with those of Zhurbas et al. (2014) but in terms of a different method.

## 3. The PDF of residual velocity

According to Taylor’s (1922) classical theory, the diffusivity can be derived from the Lagrangian integral time scale along with the mean square of velocity fluctuation, assuming that the eddy field is statistically homogeneous and stationary. A less stringent condition was assumed in Davis’s (1987, 1991) theory, that is, the probability density function (PDF) of residual velocity should be approximately Gaussian in distribution. Therefore, it is necessary to investigate whether the PDF of residual velocity observed by Lagrangian drifters is approximately Gaussian or not before exploring the diffusivity as well as Lagrangian time and length scales in the IO.

Residual velocity defined here is different from previous studies. Bracco et al. (2000) first grouped drifter observations into bins and then define the residual velocity with respect to the bin mean. Removing the mean of each bin can greatly increase the degree of freedom (LaCasce 2008) and thus the accuracy of the PDF of residual velocity. In the present study, residual velocity is defined as the departure from a mean that is estimated by the GM method. The GM mean could be a time-independent constant within a bin or include the annual/semiannual cycle, yielding different residual velocities. The main difference is whether to take the annual/semiannual cycle as a part of the mean component. A final step of dividing by the standard deviation of each bin is applied to obtain the standardized residual velocities, which is similar to Bracco et al. (2000).

Figure 1 shows the PDFs of the residual velocities after excluding the mean and the seasonal cycles over the IO, as well as the fitted Gaussian distribution (solid curve) and the double exponential distribution (dashed curve) that are calculated following Bracco et al. (2000):

where the parameters and are the standard deviation and mean obtained by a maximum likelihood estimate of the zonal or meridional residual velocity component. In general, the PDFs of both zonal (left column of Fig. 1) and meridional (right column of Fig. 1) residual velocities agree with the Gaussian distribution more than the exponential distribution, deviating only slightly (significantly) from the Gaussian (exponential) distribution at larger residual velocities. Notice that the PDF of original velocities (after wind-slip correction) is close to exponential distribution (PQL1), indicating that removing the mean and seasonal cycles of velocity would make the velocity PDF more Gaussian. The Kolmogorov–Smirnov (K–S) test (Press et al. 1992) is employed here to find out whether these PDFs are strictly Gaussian and whether the degree of freedom used for the K–S test is computed as half of the total length of the samples divided by the typical time scale (Riser and Rossby 1983) that is taken as a constant of 5 days here. Two statistics, kurtosis and skewness, are also computed for reference. It is interesting to investigate the effect of the annual/semiannual cycles on these statistics. Table 1 gives the kurtosis (the fourth-order moments), the skewness, and the K–S probability for all PDFs of the residual velocities with low-frequency components (i.e., mean, annual, and semiannual components) being removed one by one. It is found that the kurtoses computed for all PDFs are larger than that of a Gaussian distribution (i.e., a value of 3), and all skewnesses are larger than zero, implying a positive bias of all PDFs against the Gaussian distribution. It is noticeable that, in general, as more of the low-frequency components are removed, the values of the kurtosis and skewness (K–S probability) become smaller (larger) (except the zonal K–S probability after removing all seasonal cycles), which means that the low-frequency components in the residual velocities lead to significant deviation of their PDFs from the Gaussian distribution. Such small K–S probabilities (Table 1) indicate that the PDFs of the residual velocities are different from Gaussian. However, given that their kurtoses (3 ~ 4) are not far from that (~3) of a Gaussian distribution, the PDFs of residual velocities can still be considered as a quasi-Gaussian distribution (LaCasce 2008). Therefore, Davis’s diffusive theory of the advective–diffusive formalism (Davis 1987, 1991) is applicable to the IO under the less stringent condition of approximately Gaussian-distributed PDFs of residual velocities.

## 4. Diffusivity and integral scales

### a. Asymptotic behaviors for different subregions

Because the spatial distribution of Lagrangian diffusivities is highly inhomogeneous, it is of interest to investigate the Lagrangian characteristics in different subregions associated with different circulation or water masses in the IO (Fig. 2a). Considering the density of drifter observations and the characteristics of the associated flow regimes (e.g., from the EKE analyzed in PQL1), we focus on nine subregions (see the black boxes in Fig. 2b), denoted as D1 (the southeastern Arabian Sea), D2 (the western Bay of Bengal), D3 (east of Somalia), D4 (the central-eastern equatorial IO), D5 (south of the central tropical IO), D6 (the Mozambique Channel), D7 (south of Madagascar), D8 (the south subtropical IO), and D9 (west of Australia), where the Lagrangian observations are relatively dense and the flow regime is assumed to be homogeneous inside each subregion.

The zonal and meridional diffusivities estimated as a function of time lag using the pseudotrack method are shown in Figs. 3–4 for the different subregions. Here, a maximum time lag of 60 days is chosen to see their asymptotic behaviors, which is much longer than the Lagrangian integral time scale of the order of a few days. The influence of seasonal variability on the diffusivity estimate is obvious in most regions of the IO basin, especially in the northern and equatorial IO (D1–D4) where strong annual and semiannual cycles exist because of the annually reversing monsoon winds; for these regions, the diffusivities are largely overestimated if a seasonal cycle is not removed from the residual velocity (e.g., using the traditional binning method). In particular, the semiannual cycle has the largest impact on the zonal diffusivity estimate for D4 that is located on the equator, but has little effect on the meridional diffusivity estimate as the semiannual signal is not significant in the meridional velocity (see Fig. 11d in PQL1). For regions where seasonal signals are weak or absent such as D7–D9 (Figs. 3–4), removing the seasonal cycles yields very similar results to those including seasonal cycles. For most of the regions, although the diffusivities appear to grow with time in a quadratic manner (ballistic behavior) at smaller time lags (<5 days), they tend to approach a steady value (asymptotic behavior) from ~10 days, after removing the annual/semiannual cycles. These results validate the conclusion of Qian et al. (2014) that in terms of real drifters, one has to remove the seasonal cycles to obtain the asymptotic diffusivity.

It is of interest to notice the overshooting at initial time lags and the meandering behavior at later time lags in the meridional diffusivity for D4 (Fig. 4). This could be attributed to the strong zonal equatorial jets (EJs) in the equatorial IO, which advect drifters around eddies and meanders, resulting in back and forth meandering motions and thus a negative lobe in the velocity autocorrelation (Klocker et al. 2012a,b). Studies have found that a mean flow significantly suppresses the cross-stream diffusivity (e.g., Ferrari and Nikurashin 2010; Klocker et al. 2012a,b; Naveira Garabato et al. 2011); that is why the estimated is much smaller than in D4 due to the suppression of cross-stream eddy mixing by the EJs.

It is worth noting that, for D1 and D3, both the zonal and meridional components of diffusivity asymptote to constant values at relatively short time lags (~10 days). This result provides solid evidence relevant to the current controversy regarding the particle-based estimate of diffusivities with a relatively short Lagrangian velocity time series versus a tracer-based estimate (Klocker et al. 2012a,b; Riha and Eden 2011; Sundermeyer and Price 1998). Klocker et al. (2012b) concluded that particle- and tracer-derived estimates of diffusivities are equivalent, and the discrepancy between particle- and tracer-derived estimates claimed in some studies (e.g., Sundermeyer and Price 1998) mainly lies with computation of the autocorrelations that generally ignores the negative lobe in the Lagrangian velocity autocorrelation function and thus leads to overestimated diffusivities. Klocker et al. (2012a,b) claimed that the discrepancies can only be reconciled when a sufficiently long Lagrangian time series is used for the particle-based estimate (they suggest using time lags of the order of 1 month); with too short a lag, the resulting drifter-derived diffusivities can be much too large. Their results were based on a stochastic model and synthetic float data with geostrophic velocities derived from satellite-measured surface height anomalies in the South Pacific. Our results, based on Davis’s theory using GM decomposition for residual velocities and a real drifter dataset, suggest a different conclusion: the drifter-derived diffusivities can asymptote to constant values in relatively short time lags (~10 days) in some regions such as D1 and D3 of the IO when the seasonal cycles are properly excluded using the GM method.

To see the diffusive anisotropy more clearly, we can further rotate the diffusivity tensor by an angle into major and minor components (e.g., Oh et al. 2000; Rypina et al. 2012):

where . The quantity is the major (minor) component along which diffusivity reaches its maximum (minimum) so that diffusive anisotropy can be clearly shown through this rotation. Oh et al. (2000) showed that such a rotation can effectively transfer the shear dispersion into the major component, leaving the minor component as pure asymptotic diffusivity.

The values of and after removing the seasonal cycles, as well as their corresponding and , are shown in Fig. 5. For most of the subregions (D1–D2, D4–D5, and D7–D8, left and middle columns of Fig. 5), the zonal (meridional) component is quite close to the major (minor) diffusivity component. In these subregions, the preferential direction (anisotropy) of diffusion could be seen directly from zonal and meridional diffusivities, and the meridional diffusivity could be roughly viewed as the true asymptotic diffusivity except for D5 (see Fig. 4). For the subregion along the Somalia jet (D3) where the mean flow is southwest–northeast oriented, anisotropy cannot directly be seen by comparing the zonal and meridional diffusivities because they are not significantly different from each other and could be highly correlated with each other with comparable cross term , that is, . After rotating into major and minor components, the difference is enlarged and the magnitude of the major component is almost twice as large as that of the minor one, indicating an obvious diffusive anisotropy in D3. The situation in D9 is very similar to that in D3 (Fig. 5). The obvious anisotropy of diffusivity over the EJs, SC, and ARC (parts of ACC) is presumably caused by the diffusion suppression crossing mean flows or boundaries. The largest anisotropy in the EJs (D4), however, could be caused by some other factors besides diffusion suppression, for example, the equatorial waveguide that facilitates zonal dispersion.

It is noticeable that even the minor components in D5 and D6 do not asymptote to a constant; instead they increase slightly with time lag. This may be caused by nonlocal estimates for larger time lags when drifters sample different diffusion regimes as they travel farther from the domain of interest. As shown in Fig. 6, the drifter tracks used to estimate the 60-day lagged diffusivities for both D5 and D6 spread much larger than the coverage of both domains. Some drifters in D5 even travel to the equator where the zonal diffusivity would be much larger than that in D5. Drifters in D6 could also travel outside the region after 60 days. Thus, the gradual changing of the diffusion regime as drifters travel and sample could be the primary reason why diffusivities in some subregions such as D5 and D6 do not approach constants (different diffusion regimes around D6 can be seen more clearly later in Fig. 7, which shows the diffusivity ellipse at 2° resolution).

To obtain an asymptotic estimate of diffusivity, some previous studies (Andersson et al. 2011; Klocker et al. 2012b; Koszalka and LaCasce 2010) have pointed out that it is better to take an average of diffusivities over the “intermediate” period instead of a maximum, that is, choosing an intermediate averaging period after the initial transient but before the measurement errors have grown too large or the diffusion regime changes significantly. Moreover, the choice of an intermediate period should also take into account the fact that the mean flow shear may make the diffusivity estimate increase with time lag, that is, to make the estimate as little contaminated by shear as possible. Based on the asymptotic behaviors of diffusivity for most subregions of the IO shown in Fig. 5, we choose 10–15 days as the intermediate period for obtaining the asymptotic values; such an intermediate period is right after the initial ballistic increase of diffusivity when it then becomes approximately flat and thus can avoid both larger measurement errors and nonlocal estimates as time grows. On the other hand, however, this intermediate period (10–15 days) is about twice the Lagrangian integral time scales and may underestimate the diffusivities (e.g., Chiswell 2013). To assess such bias, following Chiswell (2013), we estimate the percentage of the true diffusivity value Eq. (1) returns when the 10- ~ 15-day average is used, assuming the behavior for Eq. (1) can be crudely modeled as . In the selected nine subregions, Eq. (1) returns a percentage of 96% (99%) of the true zonal (meridional) diffusivity on average, with the lowest percentage of 84% (97%) in D4 (D5) where the Lagrangian time scale is the largest (see Table 2).

Table 2 gives the Lagrangian statistics for all nine subregions, including diffusivity, time, and length scales, which are calculated using residual velocities after removing the mean, annual, and semiannual cycles. In general, the diffusivities, integral time scales, and length scales for different subregions vary in a broad range of 6.0–85 cm^{2} s^{−1}, 2–7 days, and 33–223 km, respectively. The most noticeable feature of the Lagrangian statistics is seen in the equatorial IO (D4) associated with the EJs, which has the largest (or ) (~85 × 10^{7} cm^{2} s^{−1}), (~7 days), and (~223 km) and shows the most significant anisotropy with (or ) an order of magnitude larger than (or ; also see the corresponding time and length scales). The zonal and meridional estimates in the equatorial IO (84.9 × 10^{7} and 7.6 × 10^{7} cm^{2} s^{−1}, respectively, in D4) are comparable to those in the equatorial Pacific (73.5 × 10^{7} cm^{2} s^{−1} and 9.1 × 10^{7} cm^{2} s^{−1}, respectively) estimated by Bauer et al. (2002). It is also worth noting that the magnitudes of zonal and meridional diffusivities in the eastern coast of Somalia (D3) are both large (~31 × 10^{7} cm^{2} s^{−1}) and comparable to each other with a near 45° inclination angle. The south subtropical IO (D8) where the subtropical gyre is located with less active mode water, however, has nearly the smallest values of diffusivity (~6–7 × 10^{7} cm^{2} s^{−1}), time scale (2–3 days), and length scale (~40 km).

It is worth noting that the minor principal component diffusivity is responsible for the cross-stream growth of tracer dispersion. The major principal component , approximately along stream and influenced by the coupled effect of cross-stream diffusion and the mean current shear, grows without limit and therefore eventually become larger than the “pure” along-stream diffusivity (e.g., Csanady 1973). That is why knowledge of is regarded as more important than that of with respect to the intrinsic diffusion (e.g., Oh et al. 2000).

### b. Spatial distribution in the entire IO basin

Similar to that of Zhurbas et al. (2014), a 5° × 5° bin, within which drifter observations are treated as the origins of pseudotracks, is used to estimate the asymptotic diffusivities over the entire IO. Such a bin moves with a 2° offset in both zonal and meridional components, yielding a map with 2° resolution. To avoid the initial transient phase at smaller lags and the growing errors at larger lags (Andersson et al. 2011; Davis 1991; Klocker et al. 2012b; Koszalka and LaCasce 2010), the 10- ~ 15-day mean used previously is also adopted here as the asymptotic value. Diffusivities computed by removing the GM mean, GM mean and annual cycle, and GM mean, annual, and semiannual cycles are plotted in Fig. 7 as ellipses following Rypina et al. (2012).

Generally, in coastal regions diffusivity ellipses are oriented along the coastline, showing significant anisotropy because of the inhibition of cross-coastline diffusion. This is most obvious in the Somalia Current (SC) where the mean flow is also along the Somalia coast. In the interior of the IO, significant anisotropy is found around the equatorial region. Diffusivity ellipses are all oriented along the zonal direction. This is because of the equatorial waveguide facilitating the zonal propagation of tropical instability waves (Gill 1982) and thus facilitating zonal dispersion, while the EJs along the equator tend to suppress the meridional dispersion (Klocker et al. 2012a; Ferrari and Nikurashin 2010). It is worth noting that before excluding the seasonal cycles (Fig. 7a), the diffusivity ellipses in the EJs and SC regions are much larger, with the magnitudes of major component in both regions exceeding 200 × 10^{7} cm^{2} s^{−1}. These values are comparable to those of Chiswell (2013) who estimated zonal/meridional diffusivities over the Indian and Pacific Oceans based on the half growth rate of dispersion without explicitly taking into account the seasonal variability (see their Fig. 6). By removing the annual/semiannual cycles one by one (Figs. 7b,c), the magnitude of diffusivity is greatly reduced over the EJ and SC regions, but the anisotropy is still significant. By comparison, the distributions in Fig. 7c are generally in agreement with those of Zhurbas et al. (2014). However, the major component along the equator (~100 × 10^{7} cm^{2} s^{−1}) is significantly larger than that (~16 × 10^{7} cm^{2} s^{−1}) of Zhurbas et al. (2014), and the maxima centered along the equator are not clearly seen in their results. This difference could be partially attributed to their estimating method as expressed by Eq. (4), as mentioned in section 2. The diffusivity magnitude over other regions are more consistent with Zhurbas et al. (2014), including the minor component .

It is also interesting to compare the variance ellipses (Fig. 7d) with diffusivity ellipses (Fig. 7c). Generally, they agree well with each other in patterns. One of the most obvious differences is that the variance ellipses over the Great Whirl (GW) region show some circling feature similar to the GW circulation, while the diffusivity ellipses do not. Another obvious difference is that the variance ellipses over the ARC are oriented perpendicular to the ARC while the diffusivity ellipses are parallel to the ARC. According to the equation , large diffusivity with relatively smaller variance in one direction may result in a much larger time scale that is seen in the ARC area (Figs. 7c, 8a,b). Rypina et al. (2012) thus concluded that diffusivity cannot be directly inferred from EKE over such regions because of their different anisotropic directions.

Besides the diffusivity maps, we also plot the maps of the Lagrangian integral time and length scales (Fig. 8). There is a general tendency for the zonal component of both the Lagrangian length and time scales to increase toward low latitudes, that is, () changes from 30–50 km (2–4 days) at midlatitudes to over 200 km (6–8 days) near the equator. The time scales generally agree with those of Rupolo (2007) but also show some obvious differences, mainly because Rupolo’s estimates were based on the semisum of zonal and meridional values. In several previous studies, the Lagrangian zonal length scale *L* was found to be closely related to the first-mode baroclinic Rossby radius of deformation *R*_{1} (Fig. 9). In particular, Oh et al. (2000) obtained a relationship *L **R*_{1} in the northwest Pacific and the East Sea. Zhurbas and Oh (2003), as well as Zhurbas et al. (2014), also found such a relationship held for the midlatitude region of the Pacific. To investigate this, we compute the major and minor length scales defined as and , respectively, and then compare them to *R*_{1} in the IO.

Figure 10 shows the ratios of these length scales to *R*_{1}. Similar to Fig. 6 in Zhurbas et al. (2014), the ratio of the minor length scale to *R*_{1} (Fig. 10b) is close to zero over the equatorial band but remains roughly one outside the equatorial region, indicating equality between the minor length scale and *R*_{1} at midlatitudes. This is likely because eddies at midlatitudes are dominated by the first mode of baroclinic instability, whereas eddies in the equatorial region are influenced strongly by other factors besides the first mode of baroclinic instability, such as the mean flow shear. The ratio of the major length scale to *R*_{1} (Fig. 10a) is much larger than one south of 20°S except in the low EKE region between 25° and 35°S. A significant difference in our study from Zhurbas et al. (2014) is that the major length scale in our study is roughly the same as *R*_{1} over the equatorial band (Fig. 10a). This is because our (~100 × 10^{7} cm^{2} s^{−1}) is much larger than theirs (~16 × 10^{7} cm^{2} s^{−1}) in the equatorial region. This difference should be reconciled by further studies.

## 5. Lagrangian spectra

As pointed out by Lumpkin and Flament (2001), the integral scales are essentially the first-moment estimates of the Lagrangian eddy velocities and do not give detailed information of the cyclonic or anticyclonic distribution of energy or higher-moment velocity spectra. Therefore, it is worth examining the mean Lagrangian spectra that can give such information. Similar to the method used in Lumpkin and Flament (2001), we divide all the drifter trajectories into 120-day segments with 50% overlapping and apply a 10% Tukey window to obtain the spectra. To obtain a sufficient sampling number for the mean spectra calculation while considering the potential difference of the Lagrangian spectra in different regions of the IO, we focus on the five subregions: the Arabian Sea, BoB, tropical IO, south tropical IO, and the south subtropical IO, as denoted by the gray boxes P1–P5 in Fig. 2. These regions are much larger than those used for diffusivity estimates because we want to demonstrate low-frequency spectra up to 120 days (~4 months). From Fig. 6 it is shown that drifter tracks may spread much wider after 60 days. Thus, a larger region could increase the samples of track segments and ensure statistical significance for 120-day segments.

Figure 11 shows the variance-preserving rotary spectra of the zonal and meridional components of drifter velocities, in which the area under the curve represents the variance (or equivalently energy) contribution of each frequency band to the total variance of the series (Emery and Thomson 2001). Standard error bars, which are obtained via bootstrapping, are plotted for each mean spectral curve. Note that the direction of cyclonic or anticyclonic vortices in the Northern Hemisphere is opposite to that in the Southern Hemisphere. It is found that the BoB (P2) and tropical IO (P3) are more energetic than other regions, with the most energetic eddy-containing band at ~20–30 days. In the BoB, at the dominant band (~20–30 days), the anticyclonic energy spectrum is significantly larger than the cyclonic one; at lower frequencies (>40 days), the cyclonic spectrum is slightly larger than the anticyclonic one, in agreement with Cheng et al. (2013). With a weaker energy spectrum compared to that of the BoB, the Arabian Sea shows no obvious difference between cyclonic and anticyclonic spectra. Studies have found that the equatorial-originated Kelvin waves propagate eastward to the west coast of Sumatra and then become coastal Kelvin waves that travel along the coastal rim of the BoB; the coastal Kelvin waves radiate Rossby waves at the eastern rim of the BoB that propagate westward to the central-western BoB (Sreenivas et al. 2012; Cheng et al. 2013), leading to intensive eddy activity in the central-western BoB. The dominating anticyclonic eddies are probably a result of stronger and longer-persistence downwelling Kelvin waves associated with stronger westerlies compared to the upwelling Kelvin waves. Although the coastal Kelvin waves originated in the equatorial IO can travel to the eastern coast of Arabia Sea, they are much weaker than those around the rim of the BoB, resulting in less eddy activity in the Arabian Sea compared to that in the BoB. High-frequency energetic bands (1–3 days) are found in the off-equatorial subregions (i.e., P1, P2, P4, and P5), which are anticyclonic in both hemispheres and associated with near-inertial oscillations. Both the south tropical IO (P4) and south subtropical IO (P5) show lower energy spectra, but with a higher-frequency shift of the energetic eddy-containing band from P4 to P5, as the period of the near-inertial oscillations decreases with latitude from ~2 days in P4 to ~1 day in P5; these results are in agreement with the spatial distribution of diffusivities and integral time scales (see Figs. 7–8). The standard errors in P4–P5 are much smaller than those in P1–P3, primarily because the number of drifter observations in P4–P5 is much larger than that in P1–P3.

To see the spectral slope that indicates the energy cascade between different frequency bands, we plot the power spectra of surface velocities in log–log format for the five subregions (Fig. 12). In P1 and P2, a slope of −2 at the intermediate frequency band is found, similar to that found by de Verdiere (1983) and Lumpkin and Flament (2001), who attributed it to an off-resonant response to direct white-noise wind forcing. A steeper slope of −2.5 is found in P3, while a shallower slope of −1.5 occurs in P4 and P5. According to Rupolo et al. (1996), we have the following relationship

where and are the structure functions of low- and high-frequency bands, and *n* denotes the absolute value of the Lagrangian spectral slope. Therefore, a steeper slope (*n* > 2) implies that high-frequency motions have a small impact on dispersion (as in P3), whereas a shallower slope (*n* < 2) suggests high-frequency motions have a large impact on dispersion (as in P4 and P5).

## 6. Summary

Lagrangian statistics in the Indian Ocean are investigated using the GDP drifter dataset updated through June 2013. To obtain more realistic Lagrangian statistics, the spatial shear and annual/semiannual variability are removed in the calculation of residual velocities using the Gauss–Markov (GM) decomposition (Lumpkin and Johnson 2013). Then the Lagrangian diffusivities and integral scales are computed following the method of Davis (1991) and Poulain (2001).

The probability density functions of the residual velocities are found to follow a quasi-Gaussian distribution, which satisfies the less stringent condition required by Davis’s diffusive theory of the advective–diffusive formalism. Therefore, Davis’s diffusive theory is applicable to the IO.

The Lagrangian statistics are spatially inhomogeneous and anisotropic over the IO basin and are significantly affected by the annual/semiannual variability. The annual and semiannual variability have large impacts on the asymptotic behaviors of diffusivity estimates in the northern IO and equatorial IO, respectively, while the impacts are negligible south of 10°S. After removing the seasonal variability, asymptotic behaviors are easily found over all regions except the SEC and Mozambique Channel. The nonasymptotic feature of diffusivity for the SEC and Mozambique Channel regions could be attributed to the inhomogeneous diffusion regimes sampled by the drifters as they move. To avoid the initial transient phase at smaller lags and the growing errors at larger lags, as well as to minimize the estimated errors due to the changing of the diffusion regime, we choose averages of 10–15 days as the final estimates of the diffusivity components. The diffusivities, integral time scales, and length scales in different subregions of the IO vary in a broad range of ~6–85 × 10^{7} cm^{2} s^{−1}, ~2–7 days, and ~33–223 km, respectively. Large diffusivities (>20 × 10^{7} cm^{2} s^{−1}) occur in the equatorial IO and the eastern coast of Africa, whereas small diffusivities (~6–8 × 10^{7} cm^{2} s^{−1}) appear in the subtropical gyre of the southern IO and the north/east portions of the Arabian Sea. In particular, the equatorial IO associated with EJs has the largest zonal diffusivity (~85 × 10^{7} cm^{2} s^{−1}), corresponding to the largest time scale (~7 days) and length scale (~223 km). It is also found that the minor length scale is approximately equal to the first-mode baroclinic Rossby radius of deformation *R*_{1} at midlatitudes (*R*_{1} ~ 30–50 km), while at low latitudes the major length scale is about equal to *R*_{1}*.*

With respect to the current controversy regarding drifter-based diffusivity estimate with relatively short Lagrangian velocity time series versus tracer-based estimate, the present study reports a significant result: the drifter-derived diffusivities indeed asymptote to constant values in relatively short time lags (~10 days) for some of the IO when the seasonal cycles are properly excluded using the GM method, compared to the longer time lags of the order 1 month that are claimed to be required in some previous studies. This is an important contribution to the ongoing debate.

The Lagrangian rotary spectra indicate that anticyclonic motions dominate in most subregions of the IO basin. In the off-equator subregions, there are two energetic eddy-containing bands, that is, one-to-several weeks and one-to-several days, with the latter associated with near-inertial oscillations; in the equatorial IO, the strong, energetic eddy-containing band is from one week to several months, peaking at about two weeks. In particular, the Bay of Bengal is dominated by anticyclonic eddy motions with a period of ~20–30 days.

It is worth noting that interannual variability, such as the ENSO and Indian Ocean dipole (IOD), may have an influence on our estimates of Lagrangian statistics in the IO. The interannual variability, however, cannot be easily diagnosed by surface drifters alone because of the large variation of the number of observations from one year to the next. Therefore, the present study cannot analyze the influence of interannual variability on Lagrangian statistics. Finally, as pointed out by Lumpkin and Flament (2001), caution should be taken when applying the estimated diffusivities in a numerical model due to the potential presence of a red cascade of eddy energy such as in regions of eddy to mean momentum convergence that may lead to a negative eddy diffusivity.

## Acknowledgments

The authors thank all operational agencies and researchers who deployed drifters in the Indian Ocean and made this study possible. This work was jointly supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA11010304), the MOST of China (2011CB403505 and 2010CB950302), National Natural Science Foundation of China (41376021 and 41306013), the Knowledge Innovation Program of the Chinese Academy of Sciences (SQ201305), Chinese Academy of Sciences through the project KZCX2-EW-208, and the Hundred Talent Program of the Chinese Academy of Sciences. R. Lumpkin was funded by NOAA’s Climate Program Office and the Atlantic Oceanographic and Meteorological Laboratory.

## REFERENCES

^{3}He

*Numerical Recipes in Fortran 77: The Art of Scientific Computing.*2nd ed. Cambridge University Press, 933 pp.