Stirring in the Tasman Sea is examined using surface geostrophic currents derived from satellite altimeter measurements. Calculations of the distribution of finite-time Lyapunov exponents (FTLEs) indicate that the stirring in this region is not uniform and stretching rates over 15 days vary from less than 0.02 day−1 to over 0.3 day−1. These variations occur at both small (∼10 km) and large (∼1000 km) scales and in both cases are linked to dynamical features of the flow. The small-scale variations are related to the characteristics of coherent vortex structures, and there are low FTLEs inside vortices and filaments of high FTLEs in strain-dominated regions surrounding these vortices. Regional variations in the stirring are closely related to variations in mesoscale activity and eddy kinetic energy (EKE). High values of mean FTLE occur in regions of high EKE (highest mean values of around 0.2 day−1 occur in the East Australia Current separation region) whereas small values occur in regions with low EKE (mean values around 0.03 day−1 in the east Tasman Sea). There is a compact relationship between the mean FTLEs and EKE, raising the possibility of using the easily calculated EKE to estimate the stirring. This possibility is even more intriguing because the FTLE distributions can be approximated, for the time scales considered here, by Weibull distributions with shape parameter equal to 1.6, which can be defined from the mean value alone.
Understanding horizontal dispersal in the oceans is important for a wide range of problems, including plankton dynamics (Martin 2003), larval transport (Bradbury and Snelgrove 2001; Chiswell et al. 2003), the fate of pollutants (e.g., Hazell and England 2003), spatial distribution of temperature and chlorophyll (Abraham and Bowen 2002; Mahadevan and Campbell 2002), and development of eddy parameterizations (Müller and Garrett 2002). Qualitatively, it is known that horizontal dispersal of a tracer patch occurs through a combination of advective stirring, stretching and folding a coherent patch into filaments, and diffusion, mixing the patch with outside water (Eckart 1948; Garrett 1983). This stirring of tracers and production of distorted filaments is clearly seen in the evolution of tracers from deliberate release experiments (e.g., Abraham et al. 2000; Ledwell et al. 1993) and also in satellite imagery of sea surface temperatures and ocean color (e.g., Fig. 6 of Müller and Garrett 2002). However, the spatial and temporal variations in the stirring and mixing and the relationships with the flow dynamics are not well understood or quantified. Here we use velocity fields derived from satellite altimeter data to examine the spatial and temporal variations in the horizontal stirring in the surface ocean, and the relationships between this stirring and dynamical features of the flow.
Information on the horizontal dispersion in the surface ocean can be obtained from calculations of lateral diffusivities from the dispersion of surface drifters (e.g., Lumpkin et al. 2002; Zhurbas and Oh 2004, and references therein). However, the characterization of dispersion in terms of standard diffusion coefficients is appropriate only in the long-time (asymptotic) limit when particle separations are much larger than the characteristic length scale of the velocity field and particles feel uncorrelated velocity fields. This is not typically the case in the surface ocean for time scales less than two or three months, and these drifter-derived diffusivities do not accurately characterize the stirring process and production of filamentary structures. To quantify the stirring and dispersal over weeks to several months, nonasymptote diagnostics need to be used (Artale et al. 1997; Lacorata et al. 2001). Such diagnostics include the mean square separation of particles (“relative dispersion”), the mean square distance of particles from initial positions (“absolute dispersion”), the finite-time Lyapunov exponent (FTLE), and finite-size Lyapunov exponent (FSLE). [See Boffeta et al. (2001) for comparisons between these diagnostics.]
Several recent studies have used these diagnostics to quantify the surface dispersion in different ocean regions. For example, Abraham and Bowen (2002) calculated FTLEs for the East Australian Current separation region using surface currents from a blending of sea surface velocity derived from sequential satellite images of sea surface temperatures and lower-resolution data derived from satellite altimeter measurements of sea surface height (Wilkin et al. 2002), Lacorata et al. (2001) calculated relative dispersion and FSLEs for the Adriatic Sea using drifter trajectories, LaCasce and Ohlmann (2003) calculated these diagnostics plus FTLEs using drifters in the Gulf of Mexico, and d'Ovidio et al. (2004) calculated FSLEs for the Mediterranean Sea using modeled currents. Consistent with the above discussion, these studies found that over the space and time scales of several hundred kilometers and several months the dispersion was not well modeled by standard diffusion and that the mesoscale dispersion (scales less than 50–100 km) was characterized by exponential separation with e-folding time between 2 and 10 days (see also Provenzale et al. 1991; Sanderson and Booth 1991).
The above calculations of the Lagrangian dispersion have been restricted to limited regions and periods where there are suitable currents or trajectories (e.g., sufficiently high density of drifters). However, surface currents can be derived from satellite altimeter measurements of sea level (see the next section), and these measurements are available globally for over a decade. This raises the possibility of using these altimeter-derived currents to quantify the spatial and temporal variation of the Lagrangian dispersion in the surface oceans. In this case study we use altimeter-derived currents to quantify the stirring in the Tasman Sea and east of New Zealand; see Fig. 1.
Although all of the above diagnostics could be calculated using altimeter-derived currents, in this case study we calculate only finite-time Lyapunov exponent (FTLEs). FTLEs characterize the exponential stretching rate of neighboring particles and can be interpreted as the lengthening rate of a small patch of tracer (see next section). Furthermore, there exist theoretical relationships between FTLEs and the gradients of passive tracers (e.g., Pierrehumbert and Yang 1993; Neufeld et al. 2000; Abraham and Bowen 2002). This means that calculations of FTLEs can be used to understand/quantify the dispersal of (deliberate or accidental) chemical releases and to interpret structure observed in sea surface temperature and ocean color images. FTLEs have been used to quantify stirring in many different flows, including simulations of two-dimensional turbulence (Lapeyre 2002; Bracco et al. 2004) and a large range of geophysical flows, including flow in the oceans (e.g., Yang 1996; Abraham and Bowen 2002; Joseph and Swathi 1999), the atmosphere (e.g., Pierrehumbert and Yang 1993; von Hardenberg et al. 2000; Bowman 1993; Ngan and Shepherd 1999), and the mantle (Farnetani and Samuel 2003). Use of this diagnostic in our study will enable comparison of the stirring in our study region with these other flows.
The altimeter data is available globally, but we focus here on one particular region to enable a more detailed analysis. The Tasman Sea was chosen for several reasons. First, there is a complex flow pattern with a wide range of mesoscale activity within this region (e.g., Ridgway and Dunn 2003; Tilburg et al. 2001). An intense western boundary current, the East Australia Current (EAC), flows southward close to the east Australian coastline, and there is strong mesoscale activity and formation of eddies where this separates from the Australia coast around 32°S (e.g., Nilsson and Cresswell 1981). In contrast, there are very low levels of mesoscale activity in the eastern Tasman Sea and south of Chatham Rise. This wide range of mesoscale activity within the study region enables the dependence of stirring on the level of activity to be examined. Second, the stirring in this region is of oceanographic interest on its own, for example, for understanding larval transport (Chiswell et al. 2003) and ocean color (Tilburg et al. 2002). A final reason is that a preliminary examination of stirring in the EAC separation region was performed by Abraham and Bowen (2002) using a different surface current dataset. Comparison of our calculations with those of Abraham and Bowen (2002) will provide information on the sensitivity to the source of the ocean current data.
In the next section the velocity data and diagnostics used in this study are described. Then in section 3 the spatial and temporal distributions of FTLEs are described. The relationships between stirring and dynamics quantities is then examined in section 4. In section 5 we discuss the possible impact of the limited resolution of the altimeter-derived currents. Conclusions and discussion of future work are in the final section.
2. Data and methods
a. Velocity data
The surface ocean currents used in this study are derived from satellite altimeter measurements of sea level. Gridded sea level anomalies from merged Ocean Topography Experiment (TOPEX)/Poseidon and European Remote Sensing Satellite-1 and -2 (ERS-1/2) measurements were obtained from the Archiving Validation and Interpretation of Satellite Data in Oceanography (AVISO), and surface velocities were then calculated from these fields assuming geostrophy. Because of uncertainties in the geoid, the altimetric height measurements provide information only on deviations of sea level from a long-term mean, and thus the calculated velocities are also only variations from the time mean velocities. To obtain the full time-varying velocity fields we use data from the Commonwealth Scientific and Industrial Research Organisation (CSIRO) Atlas of Regional Seas (CARS; Ridgway et al. 2002) (see Fig. 1). A climatological-mean velocity is derived using the climatological-mean temperature and salinity from CARS relative to a level of no motion at 2000 m. These time-mean velocities are then added to the time-varying velocities calculated from the AVISO sea level anomalies. [Several calculations were repeated using the climatological-mean flow from Rio and Hernandez (2004) instead of CARS to form the total velocity fields, and these produced very similar results to those presented below. This suggests that there is only weak sensitivity to the mean flow used.]
In deriving the above velocities we have assumed geostrophy, and there is no account of Ekman flows or other ageostrophic currents. These affects will play some role in the stirring and dispersal of tracers, but we focus here only on the stirring due to geostrophic flows and leave examination of the impact of Ekman and ageostrophic flows to later studies.
The velocity dataset formed from the above procedure covers the 9-yr period between 1993 and 2001 with data every 7 days on a ⅓° Mercator grid. In our analysis we focus on the region between 25° and 45°S and between 150° and 185°E; see Fig. 1. The original gridded sea level anomalies are formed from the altimeter measurements using a mapping procedure with space and time correlation scales of around 150 km and 15 days, respectively (Ducet et al. 2000), so the effective resolution is lower than that of the gridded data. The possible impact of this limited resolution on our calculations is discussed in section 5 below.
b. Finite-time Lyapunov exponent
To quantify the stirring in the above velocity fields we calculate the distributions of FTLEs, defined as
where δx(τ) is the separation after time τ between two points that were initially close together [i.e., the initial separation δx(0) is small] and the orientation of δx(0) is chosen so that λ is maximal. An initially small, circular patch of tracer will be deformed into an ellipse, and λ is the rate of increase of the major axis of this ellipse. (In a pure-strain flow the contour remains elliptical, and the stretching rate λ is equal to the strain rate.)
Spatial variations in λ characterize the existence of regions with different stirring. A large value of λ corresponds to a location where there is rapid, exponential separation of two nearby particles (so-called chaotic stirring), whereas for small values of λ there is slow or no exponential separation of particles (weakly chaotic or regular stirring).
The magnitude of λ also varies with integration time τ. As τ increases there is, in general, a larger possible range of straining rates sampled by parcels and less likelihood that parcels will experience extreme values over the whole trajectory. So, in general, the range of λ within a domain decreases with increasing τ. In fact, for closed, ergodic systems the exponents converge to a single value; that is, limτ→∞λ(x) = λ∞, where λ∞ is “the” Lyapunov exponent of the flow system.
In this study we calculate the FTLEs as in Abraham and Bowen (2002). Particle trajectories are calculated using a fourth-order Runge–Kutta integration with a daily time step, with velocities interpolated in space and time to the particle locations. The FTLE along each trajectory is then calculated from the eigenvalues of the integrated deformation matrix; see Abraham and Bowen (2002) for details. (Note, Abraham and Bowen referred the FTLEs to the end point of the trajectories whereas we refer the FTLEs to the initial point.) We calculate FTLEs, to τ = 30 days, for trajectories initialized on a regular latitude–longitude grid for different dates between January 1993 and December 2001. Several series of calculations were performed for different grid resolution and domains. The first series to be discussed used a 0.5° by 0.5° latitude–longitude grid covering the whole study region, and 12 calculations (initialized at the start of each month) were performed per year for each year between 1993 and 2001. This set of calculations provides an overview of the spatial and temporal variations of the stirring in the EAC region. Later calculations use higher resolution over smaller regions (e.g., “EAC separation” region: 28°–42°S, 148°–160°E) to examine smaller-scale variations in the stirring.
c. Comparison with stretching of material contours
One possible interpretation of FTLEs is the rate of lengthening of the circumference of a small, circular material contour. If such a relationship exists, then FTLEs not only characterize regions with different stirring but also provide estimates of the stretching of patches of tracers. To explore the relationship we perform direct simulations of the evolution of material contours. Comparison of the lengthening of these contours with the FTLEs at the center of the initial contour tests the above relationship. Furthermore, the evolution of the contours illustrates the stirring in the flow.
The evolution of the material contour is simulated using the contour advection (CA) technique (Waugh and Plumb 1994; Norton 1994). In this method the contours are represented by a series of particles, which are advected by a specified flow, and the resolution of the contours is preserved by continually adjusting the number of particles (i.e., new particles are added when the separation of neighboring particles exceeds a specified maximum value). Using the CA method the evolution of small circular contours (radius ∼11 km) are simulated for 30 days, and the rate of exponential increase of the contour length calculated; that is, the length of contour l(t), is fit by l(t) = l0eSt. [Note that the rate of lengthening of material contours has been used in many studies of the stratosphere to estimate stretching rates (e.g., Pierce and Fairlie 1993; Waugh et al. 1994; Waugh and Rong 2002). However, in these calculations the initial contours were of much larger scale and circled the globe, either as zonal contours or following the polar vortex edge.]
Figure 2a shows the evolution of two contours initialized on 19 November 1997 at two different locations, marked A (36°S, 156°E) and B (35°S, 153°E). There is a dramatic difference in the evolution of the two contours: Contour A is stretched out into a long thin filament, becomes convoluted, and part of the contour folds back on itself, whereas contour B remains small and elliptical. Location A is in a region of chaotic stirring whereas B is in a region of weakly chaotic or regular stirring. The lengthening of these two contours is quantified in Fig. 2b. As expected from the maps in Fig. 2a, there is a large difference in the lengthening of the two contours: The length of contour A increases rapidly with rate of exponential increase S ≈ 0.14 day−1 (over days 5–30), whereas the length of contour B oscillates and only increases slightly over the 30 days (S ≈ 0.015 day−1). These values of S are in good agreement with the 30-day λ values of 0.15 and 0.015 day−1 at locations A and B. As will be discussed further below, the variations in the stirring in these two locations can be related to their location relative to coherent vortices: location B is inside a vortex whereas A is in a high strain location outside of a vortex.
Similar calculations have been performed at a number of other locations. Figure 3 shows a scatterplot of contour lengthening rate S versus FTLE λ for these calculations. There is some scatter, especially for larger stretching rates, but overall there is reasonable agreement between λ and S. Thus λ can be used to estimate the lengthening of a patch of fluid in a given location, and calculations of λ should provide some guidance of deformation of the water marked in deliberate tracer releases.
d. Flow diagnostics
As well as quantifying the stirring, we wish to relate the stirring to the dynamics of the flow. We therefore compare the FTLEs with several (Eulerian) dynamical diagnostics. A commonly used diagnostic of the mesoscale variability in the oceans is the eddy kinetic energy (EKE) of the flow. The EKE (per unit mass) is given by
where u′ and υ′ are the time-varying anomalies in the zonal and meridional velocities (i.e., deviations from the long-term mean) and angle brackets denote time average. Maps of EKE have been calculated from altimeter data (e.g., Stammer 1997; Ducet et al. 2000) and from drifters (e.g., Zhurbas and Oh 2004).
We also compare the FTLEs with the strain rate
and the vorticity
where (u, υ) is the total velocity field. The initial stretching of fluid elements is given by the local strain rate, and we thus expect the FTLEs to be related to the strain rate γ. We also expect some relationship between the FTLEs and the vorticity ω, as previous studies of geophysical flows and two-dimensional turbulence have related variations in the stirring and particle dispersion to coherent (vortex) structures (e.g., Elhmaïdi et al. 1993; Koh and Legras 2002; Lapeyre 2002; d'Ovidio et al. 2004).
This parameter has been used previously in studies of two-dimensional turbulence (e.g., McWilliams 1984; Brachet et al. 1988; Elhmaïdi et al. 1993), the stratosphere (e.g., Haynes 1990; Fairlie et al. 1999; Waugh and Rong 2002), and recently the Mediterranean Sea (Isern-Fontanet et al. 2003, 2004) to partition the flows into different regions. Here Q measures the relative contribution of strain and rotation on fluid elements and can be used to identify regions of the flow where strain dominates (Q > 0) and regions of the flow where rotation dominates (Q < 0). Different “mixing” is expected in each of these regions. Okubo (1970) and Weiss (1991) showed that, if the velocity gradients are slowly evolving with respect to the vorticity gradients in a Lagrangian frame, tracer gradients evolve as exp(±Q1/2t). Thus, in regions with positive values of Q fluid elements are stretched, and there is exponential divergence of nearby particles (i.e., chaotic stirring). However, inside vortices where rotation dominates and Q < 0 there is no exponential separation and only weak (or no) stirring.
The above relationships relating tracer gradients and stretching with Q are based on the assumption of slowly varying velocity gradients, and simple counterexamples exist where these results do not hold (e.g., Pierrehumbert and Yang 1993). Extensions of the Okubo–Weiss relationship have been developed that take into account the rotation of principle strain axes (e.g., Hua and Klein 1998; Lapeyre et al. 1999). We focus here on the simpler Okubo–Weiss parameter, which does not require calculation of Lagrangian acceleration terms.
3. Stretching rates
a. Spatial distributions
We first consider the spatial distributions of λ. Figure 4a shows the distribution of λ for τ = 15 days for a single calculation initialized on 19 November 1997. (In all FTLE maps shown λ is plotted according to the starting point of the trajectories.) This map shows that there is a large range of λ within the domain with λ varying from less than 0.02 to over 0.3 day−1 (which corresponds to e-folding times varying between 50 and 3 days). This variability can occur at very small scales with large differences in λ occurring for neighboring initial locations (e.g., λ at grid points separated by 0.5° can vary by over 0.3 day−1). There is, in fact, variability at much smaller scales than resolved in this figure (see Fig. 11a below). The cause of these small-scale variations in λ is examined in section 4b.
Closer examination of Fig. 4a shows that there are also variations in λ over larger (regional) scales. These regional variations can be clearly seen in the ensemble mean ( is the mean λ at each location over the series of 108 calculations with starting date varying between 1993 and 2001); see Fig. 4b. There are high along the east Australian coast ( ≈ 0.2 day−1), but low values (around 0.05 day−1) in the eastern Tasman Sea and along and south of the Chatham Rise. The signature of the dominant dynamical features in the domain can be seen in the map; for example, there are large values in the EAC separation region, and moderate values extending from along the Tasman Front to the northern tip of New Zealand. These regional variations in the stirring are examined further in section 4a.
The maps of λ show that stirring is not uniform and there is a wide range of stretching rates within the study region. To quantify the range of λ and the dependence on the integration time τ, we calculate probability distributions P(λ, τ). Figure 5a shows P(λ) for several τ (these distributions are calculated over λ from all calculations and all grid points). The distributions for small τ are asymmetric and very broad, confirming the wide range of λ shown in the maps. As the integration time τ increases, the distributions narrow and the peak moves to lower values. This is because for longer trajectories a larger range of straining rates is sampled by parcels and there is less likelihood that parcels will experience extreme values over the whole trajectory. The variation of the mean and standard deviation of λ with τ is shown in Fig. 5b. As expected from Fig. 5a, both quantities decrease with τ. This narrowing of the distribution with longer integrations means that the spatial differences discussed above are reduced as τ increases. However, clear differences between regions are still found with τ = 30 days (see Fig. 6b).
The distributions shown in Fig. 5 are similar to those calculated by Abraham and Bowen (2002) using different velocity data for a smaller region covering the EAC separation region. Abraham and Bowen performed calculations out to τ = 80 days and showed that, while the distributions continued to narrow, the convergence was very slow and the distribution had yet to converge to the Lyapunov exponent λ∞ by τ = 80 days. The narrowing of the P(λ) distributions, shifting of peak (and mean) to lower values, and slow convergence as the integration time τ increases has also been found in simulations of tropospheric flow (Pierrehumbert and Yang 1993) and two-dimensional turbulence (Lapeyre 2002).
As noted above, there can be large differences in λ between different regions of the domain considered. These differences occur not only in the mean values within the regions but also the range of values. For example, Fig. 6a compares P(λ) for the EAC and east Tasman Sea subregions (shown by the white boxes in Fig. 4). There is a much larger mean and range (standard deviation) for λ in the EAC region. Similar differences in the distributions between the regions are observed for different τ although, as discussed above, the distributions are narrower for larger τ ; see Fig. 6b. The rate of decay with τ of the mean and standard deviation of λ varies between subregions with faster (slower) decay for regions with larger (smaller) mean FTLE.
b. Temporal variations
The above analysis has shown that there are large spatial variations in the stirring in the oceans around eastern Australia and New Zealand. We now examine whether there are also temporal (seasonal or interannual) variations in this stirring.
Figure 7 shows the variation of 15-day λ between 1993 and 2001 for averages over the whole region, the EAC separation region, and the east Tasman Sea. The thin curves correspond to λ for each calculation (initialized at the first of each month) while the thick curves correspond to the annual mean λ (averages over 12 calculations per year). This plot shows that, as compared with the spatial variations, there are only weak temporal variations in λ. For all times the stirring in the EAC separation region is much stronger than the mean value over the study region, and the stirring in the east Tasman Sea is weaker than average. The interannual variations are small for the means over the whole regions and over the subregions. Larger variability occurs between individual calculations, especially for the separation EAC region. However, there is no systematic seasonal variation, and the variations between the 9-yr means for each month are also very small (not shown).
Examination of probability distributions, for the whole and subregions, also shows only weak temporal variations in the distributions, and the above results for ensemble-mean distributions also hold for individual calculations.
c. Approximate distributions
The λ distributions shown above (as well as distributions of strain rates γ) are well fit by Weibull distributions. Weibull distributions are commonly used in investigations of wind speed distributions (e.g., Conradsen et al. 1984), and are defined by
where a and b are two free parameters. The shape parameter b of the Weibull distributions that best fit the λ distributions varies from around 1.5 for τ = 5 days to around 1.65 for τ = 30 days, but all λ distributions can be well modeled with b = 1.6; see thick solid curves in Fig. 8. The mean of the Weibull distribution is given by μ = aΓ(1 + 1/b), where Γ is the gamma distribution (for b = 1.6, we have μ = 0.9a). Hence, for fixed b the parameter a can be determined from the mean of the distribution. This has important implications as it means that, given the mean λ over a region, the distribution of λ can be estimated [using Eq. (6) with a = μ/0.9 and b = 1.6], and the differences in λ distributions between different regions can be characterized by differences in the mean λ alone.
Some insight into why Weibull distributions fit the λ distributions can be gained by considering the special case of b = 2, which corresponds to a Rayleigh distribution. Although the fit is not as good as for b = 1.6, a Rayleigh distribution is still a reasonable approximation to the λ distributions; see dashed curves in Fig. 8. [For a Rayleigh distribution,
and the peak of the distribution occurs at 2−1/2a.] A Rayleigh distribution is the probability distribution function of the magnitude of a two-dimensional vector whose two components are Gaussian distributed with zero mean and the same standard deviation 2−1/2a. Hence, if the components of the strain vector satisfied these conditions, the strain rate distribution would be a Rayleigh distribution and, from the above arguments relating the FTLEs to strain rates, so would the λ distributions. The fact that the strain rate (and λ) distributions have b slightly less than 2 indicates that these conditions are not exactly met. Analysis of the distributions of the components of the strain vector shows that the two components have very similar, symmetrical distributions with means close to zero, but their shape is non-Gaussian (the distributions decay more rapidly than Gaussian distributions). This deviation of the components from Gaussian explains the deviation of the strain rate (and λ) distributions from a Rayleigh distribution.
Although Weibull distributions with b = 1.6 are good fits to the FTLE distributions for the range of τ considered here, we expect the value of b to change for much longer integration times. In particular, for very long times it is expected that the distribution P(λ) will become Gaussian (e.g., Abraham and Bowen 2002). As a Weibull distributions with b = 3.6 is close to a Gaussian distribution, we then expect b → 3.6 as τ → ∞.
4. Relationship with dynamical quantities
In the above sections we have shown that there are large spatial, and much smaller temporal, variations in the stirring in east Australian–New Zealand region. We now examine the relationship between the stirring and the flow dynamics and what aspects of the underlying flow produce the variations in the stirring. We first consider the regional-scale variations in the time-averaged stirring and dynamics and then consider variations on smaller scales.
a. Regional-scale variations
As discussed above, the stretching rate is related to the strain rate experienced by the particles. We therefore expect regional-scale variations in the mean FTLEs to be related to the variations in the time-mean strain rate γ. Figure 9a shows that there are large spatial variations in γ in the study region. There are large values along the east coast of Australia with the maximum around 32°–35°S where the EAC separates from the coast and there is intense eddy formation (Nilsson and Cresswell 1981), and moderate values extend from the separation region to the northern tip of New Zealand (i.e., along the Tasman Front) and down the eastern coast of North Island, New Zealand. In contrast, there are low values in the eastern Tasman Sea. As expected, this map of γ bears a strong resemblance to the ensemble mean FTLE , shown in Fig. 4b, with large (small) occurring in regions with large (small) γ.
The relationship between mean strain rate and mean FTLE is quantified in Figs. 10a and 10b, which show scatterplots of 5-day and 15-day versus γ. The dots in each plot correspond to individual grid points, while the boxes correspond to averages over all grid points in 5° by 5° subregions. These plots show that the –γ relationship is reasonably compact, that is, there is small variation in for given γ. For small strain rates ≈ γ (as discussed earlier) but, as the mean strain increases, becomes less than γ. This is because particles remain in the large strain regions only briefly and during the integration period sample regions with smaller γ, and the integrated stretching rate is less than the initial large strain rate. The range of γ where ≈ γ holds decreases with increasing integration τ (cf. Figs. 10a and 10b). Again, this is related to range of strain rates sampled over the integration time: For longer τ the particles are more likely to sample regions with weaker strain rates. For τ ≥ 10 days, the –γ relationships can be reasonably well modeled by =αγ1/2, curve in Fig. 10b.
A close relationship between the stretching rates and eddy kinetic energy is also found. Figure 9b shows the time-mean EKE for 1993–2001, and this distribution resembles that of the ensemble mean FTLE (Fig. 4b) and the mean strain rate γ (Fig. 9a). In fact, as shown in Fig. 10c, there is a compact linear relationship between EKE1/2 and γ, which can be well fit by EKE1/2 = Lγ with L = 80 km. As EKE1/2 is the rms velocity and the strain rate scales with the ratio of eddy velocity to eddy length scale, L can be interpreted as an eddy length scale. As a result of this linear relationship, the relationships between and EKE1/2 are very similar to that between and γ; see Fig. 10d. Thus, to a good approximation ∝ EKE1/4 (for τ ≥ 10 days).
More rapid stirring and mixing is expected in regions with larger mesoscale activity. For example, a close relationship between the lateral diffusivity calculated from the dispersion of surface drifters and EKE has been found in previous studies (e.g., Zhurbas and Oh 2004). Indeed, the map of lateral diffusivity shown in Zhurbas and Oh (2004, their Fig. 1) resembles the map of mean FTLE shown here with very large values in the EAC separation region and high values extending along the Tasman Front and down the east coast of New Zealand. Furthermore, from scaling arguments simple relationships between eddy diffusivity and EKE or root-mean square of the height variability can be derived, and these have been used to form maps of diffusivities from altimetric measurements (e.g., Stammer 1998). However, the relationships between the diffusivity calculated from drifters and EKE do not appear to be universal and do not agree with the simple relationships derived from scaling arguments (e.g., Lumpkin et al. 2002; Zhurbas and Oh 2004).
Further analysis of the –EKE relationships is required, including determining whether the same relationships apply in other ocean basins. If the relationships are robust, then it will be possible to use the, easily calculated, EKE to provide an estimate of the mean λ (stirring) in a particular region. Furthermore, as discussed above, the variability of λ can then be estimated by using a Weibull distribution with parameters a = μ/0.9 and b = 1.6.
b. Small-scale variations and coherent vortices
Having examined the regional-scale variations in the mean stretching rates, we now examine the cause of the finescale variations in λ described in section 3a. We focus on the EAC separation region as the largest stirring occurs in this region, and there is also the largest finescale variability in this region (see Fig. 4a). To examine the details of the finescale structure we have repeated some of the λ calculations on higher-resolution grids, that is, a 0.1° by 0.1° latitude–longitude grid covering only the EAC separation region (28°–42°S, 148°–160°E).
Figure 11a shows the resulting λ distribution for 19 November 1997 (cf. with Fig. 4a). The range of values of λ is similar to that in the lower-resolution calculations but there is much finer scale variability. The λ field is characterized by a sea of very narrow filaments of high λ, which are intermingled with coherent regions with very low λ. The filaments are still likely thinner than resolved in this simulation (e.g., in many cases the differences in λ between neighboring grid points exceeds 0.3 day−1).
Previous studies, of geophysical flows and two-dimensional turbulence, have related variations in the stirring and particle dispersion to the existence and characteristics of coherent (vortex) structures (e.g., Elhmaïdi et al. 1993; Koh and Legras 2002; Lapeyre 2002; d'Ovidio et al. 2004). To examine if such relationships exist for the current flow we examine the distributions of vorticity, strain rate, and the Okubo–Weiss parameter (see section 2d). Figures 11b–d show maps of these quantities for the same date as Fig. 11a.
The vorticity field (Fig. 11b) shows numerous cyclonic and anticyclonic eddies (coherent regions of negative and positive vorticity) with the largest and strongest occurring near the Australian coast and south of 32°S. This is consistent with previous studies showing that this is an intense region of eddy formation (Nilsson and Cresswell 1981). Comparison of vorticity and strain maps (Figs. 11b,c) shows that the largest strain rates occur at the edge of these vortices with weak strain inside the vortices and in the southeastern part of the domain where there are no strong vortices.
As discussed in section 2, a measure of the relative contributions of strain and rotation is the Okubo–Weiss parameter Q, and this parameter can be used to partition the flow into strain-dominated (Q ≫ 0) and rotation-dominated (Q ≪ 0) regions. This is confirmed in Fig. 11d, which shows negative Q inside the coherent cyclones and anticyclones (cf. with vorticity field), large positive Q in regions with large strain rate surrounding these vortices, and small Q in the quiescent southeastern region.
Comparing the maps of vorticity, strain, and Q with that of λ we see that, consistent with previous studies, a strong relationship between the finescale structure in λ and the coherent vortices (the contours in Fig. 11a correspond to Q = −0.2 day−1, and identify the coherent vortices in the flow). There are, in general, very low values of λ inside the coherent vortices while the filaments of large values occur in the high strain regions surrounding the vortices.
To quantify this relationship we calculate the distributions of λ conditioned by the value of Q (Elhmaïdi et al. 1993). In particular we consider the distributions for strain-dominated regions with Q > δQ (δQ > 0) for rotation-dominated regions (vortices) with Q < −δQ and for background regions with near-zero Q (−δQ < Q < δQ). Following Pasquero et al. (2001) and Isern-Fontanet et al. (2004), we use δQ = 2σQ (where σQ = 0.1 day−1 is the standard deviation of the Q distribution in the domain). Figure 12 shows the conditional probability distribution functions (PDFs) of λ for 19 November 1997. There is a clear difference between the distributions, particularly between the distributions for strain- and rotation-dominated regions. There is a much smaller mean and range of λ in negative Q (rotation dominated) regions than in positive Q (strain dominated) regions. Similar differences between conditional distributions of λ are obtained for calculations for other dates (not shown), and the difference between stirring inside and surrounding vortices is a robust result.
Although the λ distributions in Fig. 12 clearly differ, there is still some overlap; that is, some large (small) values of λ occur with regions with negative (positive) Q. An example of high λ inside a coherent negative-Q region (vortex) can be seen (at 33°S, 156°E) in Fig. 11. The reason for this is that this vortex is sheared out and destroyed in the following 14 days (not shown). Although there is only weak strain during the initial stages of the integration, large strain occurs during the destruction of the vortex, and this results in the large λ (and stirring). Using a diagnostic that takes into account the rotation of strain axes along Lagrangian trajectories (e.g., Hua and Klein 1998; Lapeyre et al. 1999) may result in even more distinct λ distributions than shown here. However, even with the simple Q diagnostic, it is clear that the stirring within coherent vortices is much weaker than in regions surrounding these vortices.
It is interesting to compare the distributions shown in Fig. 11 with those from two-dimensional turbulence simulations. In particular, Lapeyre (2002) classified structures in the λ field calculated in a two-dimensional turbulence simulation into several types, and noted that the filamentary structures with largest exponents were around or connected between vortices, the very lowest values were inside the vortices, and there were low values in the background turbulent field. These features are also seen in our λ simulations. Thus, not only are PDFs of λ very similar to those in two-dimensional turbulence (see section 3a above), but so is the relationship between λ and coherent vortices.
Given the track spacing of the altimeters, an obvious question is: How reliable are the estimates of stirring calculated from altimeter-derived ocean currents? The answer to this is unfortunately not known. However, several pieces of information suggest that the estimates might be reasonable.
For reliable estimates of small-scale structures and stirring to be calculated from the resolvable (mesoscale) features of the flow, the larger-scale flow features must dominate the strain field. The ability of low-resolution velocity fields to produce reliable estimates of small-scale structures and stirring can be related to the steepness of the kinetic energy spectra (Shepherd et al. 2000; Bartello 2000). When there is a steep decay of the kinetic energy spectra [E(k) ∼ k−α with α > 3], there is “spectrally nonlocal” dynamics and particle dispersion is controlled by the large-scale velocity, whereas for shallow decay (α < 3) there is “local” dynamics and dispersion is controlled by the velocity at the scale of the dispersion (Bennett 1984; Babiano et al. 1985). Thus, reliable stirring estimates can be obtained from low-resolution winds having steep kinetic energy spectra but not from winds having shallow spectra. This was confirmed by the modeling studies of Bartello (2000) and Shepherd et al. (2000). In particular, Shepherd et al. (2000) showed that FTLE calculations in the stratosphere, troposphere, and mesosphere (where α is greater than, around, and less than 3, respectively) are insensitive, weakly sensitive, and strongly sensitive to the spatial resolution of the velocity field, respectively. Although there is some uncertainty in energy spectra for the surface ocean, the calculations of Stammer (1997) indicate that in the midlatitudes the energy spectra decay is around k−3. This is the boundary between nonlocal and local dynamics, and, based on the analysis of Shepherd et al. (2000), we might expect relatively weak sensitivity of FTLE calculations of the surface ocean to the spatial resolution of velocity field and reliable estimates from velocity fields with spatial resolution of the altimeter data.
Additional insight into the reliability of the λ calculations can be gained by comparison with calculations using the currents used by Abraham and Bowen (2002). The surface currents that they used were formed by the blending of velocities derived from sequential satellite images of sea surface temperatures with those from satellite altimeter measurements (Wilkin et al. 2002) and contain higher-resolution information than the altimeter-only currents used here. We have performed parallel calculations of λ for these two surface current datasets. (In these calculations the two velocity datasets were interpolated onto the same latitude–longitude grid covering the smaller domain of the blended velocities.) Although there are these point-to-point differences (which are, in general, related to differences in the shape of vortices), there is good agreement in the distributions of λ and dependence on integration time τ, and the results based on calculations using the altimeter-only currents also hold for calculations using the blended data. Although this agreement is encouraging, the resolution of the blended data is still not sufficient to resolve small-scale velocity features, and surface currents with much higher resolution than the blended velocities are needed to fully quantify the dependence of the λ calculations on the resolution of the surface currents.
Another approach to assess the reality of the calculated stretching rates is by comparison with the stretching of tracer patches from deliberate tracer release experiments. Estimates of the stirring in the surface ocean have been obtained from in situ SF6 measurements or satellite images of induced phytoplankton bloom following several surface tracer releases, including three release experiments in the Southern Ocean (Abraham et al. 2000; Coale et al. 2004) and one in the North Pacific Ocean (Law et al. 2005, manuscript submitted to Deep-Sea Res. II). The stretching rates from these experiments are between 0.07 and 0.11 day−1. A detailed comparison of these values with our estimates is not possible because of the different oceans sampled. Also, the tracer releases provide only a single estimate for each location and our calculations indicate a wide range of λ (and small-scale variability in the stretching rates). However, it is encouraging that these observational estimates are of the same magnitude as our calculations of λ. Furthermore, the stretching rates vary with the EKE at the location of the tracer release: The three releases with stretching rates around 0.07–0.08 day−1 were in a region with EKE around or less than 30 cm2 s−2, while the release with larger stretching [northern release in the Southern Ocean Iron Experiment (SOFeX) cruise] was in a location with EKE around 200 cm2 s−2.
6. Concluding remarks
Using surface currents derived from satellite altimeter measurements we have examined the stirring, as quantified by finite-time Lyapunov exponents, in the oceans east of Australia and surrounding New Zealand. These calculations indicate that the stirring in this region is not uniform, and there is a wide range of stretching rates. For example, the stretching rates over 15 days vary from less than 0.02 to over 0.3 day−1 (corresponding to e-folding times of less than 3 days to over 50 days). These variations occur at both small (submesoscale: ∼10 km) and large (regional: ∼1000 km) scales.
High-resolution maps of finite-time Lyapunov exponents, λ, show coherent regions of very low values surrounded by narrow filaments of high values. As found in other studies of geophysical flows or two-dimensional turbulence (e.g., Elhmaïdi et al. 1993; Koh and Legras 2002; Lapeyre 2002; d'Ovidio et al. 2004), these small-scale variations in λ are related to the existence and characteristics of coherent vortex structures. Comparison with the vorticity, strain rate, and Okubo–Weiss parameter shows that the low λ are in rotation-dominated regions inside coherent vortices while the filaments of high λ occur in strain-dominated regions surrounding the vortices.
Together with these small-scale variations there are also regional variations in the stirring. These variations are related to dynamical features in the region, with high values of mean λ occurring in regions of strong eddy activity (e.g., highest mean values around 0.2 day−1 occur in the EAC separation region) whereas small values occur in regions with little eddy activity (e.g., mean values around 0.03 day−1 in the east Tasman Sea). The maps of mean λ are similar to maps of lateral diffusivity derived from drifters (e.g., Zhurbas and Oh 2004) and are also very similar to maps of the mean strain rate and EKE. There is, in fact, a compact relationship between and EKE with ∝ EKE1/4. The existence of this –EKE relationship raises the possibility of using EKE, which can be easily calculated, to estimate the mean stirring within different regions. This possibility is even more intriguing as the λ distributions are well fit by Weibull distributions with parameters a = μ/0.9 and b = 1.6. Hence, given an estimate of the mean λ from the EKE in a region, it may be possible to also estimate the distributions of λ in that region.
This is a potentially important result, which may be useful for plankton and ecosystem modeling, understanding variability of passive tracers, or the development of eddy parameterizations. However, there are two important outstanding issues: 1) how reliable are estimates of stirring obtained from altimeter-derived currents and 2) do the relationships found in the Tasman Sea hold in other oceans?
The arguments and comparisons discussed in the previous section suggest that our stirring estimates may be realistic. However, explicit λ calculations using high-resolution surface currents subsampled at different resolutions (including the same sampling as the altimeters) are required to determine the exact impact of resolution on these calculations. (These calculations will also enable the impact of Ekman and ageostrophic currents on the stirring to be examined.) We plan to perform such experiments using surface currents from an eddy-permitting model.
We also plan to repeat the calculations presented here for other oceans. Preliminary global FTLE calculations have been performed, and the regional-scale variability of the FTLEs are closely related to variations in EKE in all ocean basins. However, further calculations and analyses are required to determine if the quantitative λ–EKE relationships are the same in all oceans, or whether the relationships depend on the dynamical processes controlling the flow in each ocean. It will also be interesting to examine the smaller-scale variations in the stirring and relationships with coherent vortices in the different oceans. Again qualitatively similar results are expected, but the details may vary with, for example, the vortex size and lifetime in different western boundary currents. As well as further analysis of the λ distributions, it will be of interest to perform passive tracer simulations using the altimeter-derived currents. Such calculations may provide insight into the spatial structures in satellite images of sea surface temperature and ocean color.
We thank Steve Chiswell, Mick Follows, Erik Kvaleberg, John Marshall, Graham Rickard, and Ric Williams for useful discussions and comments. The initial stages of this work were supported by a grant from the Royal Society of New Zealand ISAT fund.
Corresponding author address: Darryn W. Waugh, Dept. of Earth and Planetary Science, The Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218. Email: firstname.lastname@example.org