The robustness of Nakamura’s effective diffusivity diagnostic for quantifying eddy diffusivities of tracers within a streamwise-average framework is carefully examined in an oceanographic context, and the limits of its applicability are detailed. The near-surface Southern Ocean geostrophic flow, obtained from altimetric observations, is chosen for particular examination since it exhibits strong eddy activity. Careful analysis indicates that the effective diffusivity calculation is not strongly sensitive to the resolution of the velocity field or to the tracer initial conditions (after a spinup time of 3 months), and that a numerical diffusivity of 50 m2 s−1 is appropriate to resolve the Batchelor scale at a computational resolution of . The spatial variability of the effective diffusivity is examined and shown to be related to the eddy kinetic energy (EKE), mean zonal current (u), and phase speed of the eddies (c), following the streamwise averge of EKE/(u − c)2, in agreement with linear baroclinic instability theory. The temporal variability of the effective diffusivity over the period 1996–2001 indicates significant variability on seasonal time scales with a magnitude of ±10%–15%. Interannual variations are investigated and discussed in the context of changes in the wind stress associated with the southern annular mode.
Variability in the ocean is dominated by mesoscale fluctuations. Satellite observations indicate that baroclinic eddies are almost ubiquitous in the extratropical oceans with regions of particularly strong eddy activity being found in and around intense western boundary currents, such as the Gulf Stream and the Kuroshio, and in the Antarctic Circumpolar Current (ACC) (Stammer 1997 and references therein). These eddies transport heat and momentum globally and interact with the mean flow. They also influence the global distribution of salt, dissolved minerals/gases, and biological species. Consequently, eddy processes are crucial to determining the climate state.
Our knowledge of mesoscale eddies has substantially improved with the availability of more than a decade of global high-resolution sea surface height variability data based on satellite altimetry (Stammer et al. 2006). Analysis of these data has shed light on eddy generation mechanisms. Regions of enhanced eddy kinetic energy are seen to coincide throughout much of the extratropical ocean with the locations of mean frontal structures, suggesting the eddies extract energy from the mean flow predominantly through baroclinic instability (Stammer 1997) [this is also consistent with the findings of modeling studies: e.g., Jayne and Marotzke (2002); Best et al. (1999) for the Southern Ocean]. The eddy activity is observed to exhibit considerable temporal variability (Stammer and Wunsch 1999; Stammer et al. 2006). Meredith and Hogg (2006) hypothesize that in the Southern Ocean temporal variability in the eddy kinetic energy may be a lagged response to changes in the wind stress associated with the Southern Annular Mode.
The horizontal scale of eddies (Leddy) is observed in the extratropics to be a few times larger than the local first mode of the Rossby radius of deformation (LD), such that Leddy decreases from close to 60 km at 30° to less than 30 km at 60° (Stammer 1997; Eden et al. 2007). Thus, it would appear that the energy released by baroclinic instability from mean-flow interactions is kept close to the input wavelength;1 however, the precise details of the kinetic energy cascade for the barotropic and baroclinic modes are still debated (Vallis 2006; Scott and Wang 2005). This suggests,2 as is verified by observations, that the eddy kinetic energy in the ocean should be typically an order of magnitude larger than the mean kinetic energy.
To understand the effects of eddies on the general circulation of the ocean, it is important to understand their dynamics and quantify their mixing and transport properties. For example, the details of eddy dynamics are important for the representation of the response of the Southern Ocean overturning circulation to changes in wind stress forcing (Hallberg and Gnandesikan 2006), and near-surface diapycnal buoyancy fluxes likely play a significant role in the maintenance of the thermocline and in balancing the ocean heat budget (Marshall 1997; Marshall et al. 2002; Karsten et al. 2002; Radko and Marshall 2004; Stammer et al. 2006; Greatbatch et al. 2007).
A number of approaches have been used to quantify eddy diffusivities in the ocean, either from measurements using floats and drifters (Davis 1991; Lumpkin et al. 2002; Zhurbas and Oh 2004; Sallée et al. 2008) or by relating the eddy diffusivities to characteristics of the flow field such as eddy kinetic energy (Eden and Greatbatch 2008) or streamfunction variance (Holloway 1986; Keffer and Holloway 1988; Stammer 1998). In this paper, we will focus on quantifying eddy diffusivities and their influence on tracer distributions using Nakamura’s effective diffusivity diagnostic (Nakamura 1996). This diagnostic was used for the first time in an oceanic context by Marshall et al. (2006), who applied it to surface velocity fields derived from altimetric observations. Subsequently, Abernathey et al. (2009 hereafter AMMS) have applied the same diagnostic to full-depth model velocity fields. Both studies have provided important new information concerning the distribution of eddy diffusivities. Given the prospects for the effective diffusivity to become a standard diagnostic tool for ocean flows, it seems timely to conduct a comprehensive assessment of the robustness of the diagnostic in an oceanic context, including the sensitivity of the results to details of the calculation method and to the accuracy of the input velocity data.
Marshall et al. (2006) demonstrated the importance of the mean zonal current [not accounted for in the approaches of Eden and Greatbatch (2008) or Holloway (1986)] in determining the Nakamura diffusivity in the region of the ACC. They found much higher eddy diffusivities in this region in a study where they artificially removed the mean flow. Guided by baroclinic instability theory, here we will explore in greater detail the relationship between the Nakamura diffusivity and characteristics of the flow field by examining the causes of spatial and temporal variability.
The outline of the paper is as follows. In section 2 the diagnostic approach for quantifying eddy diffusivities is described and its robustness in an oceanographic context is analyzed using a velocity field representative of the near-surface Southern Ocean (section 3). The spatial distribution of the eddy diffusivities is examined in relation to other properties of the flow (eddy kinetic energy, mean zonal current, and phase speed of the eddies) in section 4. The temporal variability of the eddy diffusivity in the Southern Ocean is then examined in section 5 and the results are discussed in the context of changes to the wind forcing of the ocean.
2. Diagnostic approach
a. Lessons from atmospheric flows
The approach taken in this paper to quantify the eddy diffusivity of the ocean derives from atmospheric studies. In the atmosphere, eddies also play an important role in the dynamics and in transport and mixing processes (although there are some differences compared with the ocean, which are outlined below). A useful diagnostic in this context has proved to be the effective diffusivity (κeff) introduced by Nakamura (1996). Shuckburgh and Haynes (2003) demonstrated that, when applied to a variety of kinematic flows in the “chaotic advection” regime, the effective diffusivity can be used to quantify the diffusion of Lagrangian particles in a cross-streamline direction. Haynes and Shuckburgh (2000a,b), Allen and Nakamura (2001), and a number of subsequent studies (e.g., Kostrykin and Schmitz 2006) have used κeff to highlight the division in height–latitude of the upper troposphere, stratosphere, and mesosphere into eddy mixing regions separated by barriers to isentropic eddy transport in the vicinity of the stratospheric polar night jet and the subtropical jet, and have demonstrated that the effective diffusivity is a robust diagnostic in this context. However, differences described below between atmospheric and oceanic flows mean that this evidence is of only limited value in determining whether κeff can unambiguously characterize the cross-stream eddy diffusivity for ocean flows.
In the troposphere, baroclinic instability is responsible for much of the eddy activity (with the necessary available potential energy of the mean flow being replenished by radiative forcing rather than by wind and buoyancy forcing at the ocean surface); Leddy is much greater than in the ocean, being O(1000 km). This is approximately the size of LD in midlatitudes with the implication that the eddy and mean kinetic energies are comparable (unlike in the ocean) (Vallis 2006). The turbulent cascade properties display qualitative differences with respect to the ocean due to differences in the density stratification (Scott and Wang 2005). In the stratosphere, much of the eddy activity is associated with the breaking of Rossby waves of planetary scale. In contrast, in the ocean the eddy scale is much smaller than the large-scale flow; this means that there is an opportunity for eddies to be advected by mean flow into regions with quite different properties before completing their life cycle.
Many studies (e.g., Waugh and Plumb 1994; Shepherd et al. 2000) have indicated that isentropic flow in the stratosphere is in the chaotic advection regime, where it is the large-scale flow that generates small-scale features in a tracer field. We present analysis below to test whether the Southern Ocean surface velocities also lie in the chaotic advection regime and, thus, whether previous results concerning robustness of the effective diffusivity diagnostic in chaotic advection flows are relevant.
Key aspects of the robustness of the effective diffusivity diagnostic rely on the exact nature of the flow and, in particular, on the interactions between different scales of the flow. As noted above, precisely these interactions between different scales of the flow may be somewhat different in the ocean and the atmosphere. In addition, the very different time and length scales of ocean flows mean that the sensitivity of the effective diffusivity to the parameters of the calculation will be different. Here, therefore, we revisit the question of robustness of the effective diffusivity diagnostic in an ocean context. We present a careful analysis to determine the parameters within which the effective diffusivity may be considered a well-defined quantifier of the eddy diffusivity of ocean flows.
b. Effective diffusivity
Eddy fluxes in the ocean typically include a strong rotational component, which plays no role in the tracer budget. The effective diffusivity approach focuses on determining the irreversible mixing effect of eddies on tracers, which results from the divergent eddy fluxes. It does so by considering the evolution of a conservative tracer with concentration c(λ, ϕ, t) in two dimensions (λ longitude, ϕ latitude) according to an advection–diffusion equation:
where u(λ, ϕ, t) is a nondivergent velocity field and k is a constant diffusivity. The irreversible mixing effects can be isolated from the effects of advection by rewriting the equation in the form of a diffusion-only equation using a transformation to a coordinate system based on the tracer itself (Nakamura 1996; Winters and D’Asaro 1996). On a sphere of radius r, it is convenient to define the coordinate system in terms of the equivalent latitude ϕe(C, t), which is related to the area A within a tracer concentration contour c = C by the identity A = 2πr 2 (1 − sinϕe). For each tracer contour the equivalent latitude is therefore the latitude that the contour would have if it were remapped to be zonally symmetric while retaining its internal area (see Shuckburgh and Haynes 2003, Fig. 1), and the correspondence is unique, meaning C can be regarded as a function of ϕe and t.
In the new tracer-based coordinates, Eq. (1) becomes
where κeff, the effective diffusivity that quantifies the irreversible mixing, is defined by
Here Leq(ϕe, t) is the equivalent length, ϕ* is the equivalent latitude corresponding to the minimum tracer concentration contour, Lmin(ϕe) is the minimum length a tracer contour with equivalent latitude ϕe can attain given the geometry [for atmospheric flows, this is 2πrcosϕe, but for oceanic flows the presence of continental boundaries means the minimum attainable length is larger; see Marshall et al. (2006)], L represents the length of a tracer contour, and ( · ) denotes the line average around the tracer contour corresponding to ϕe. The equivalent length, and hence the effective diffusivity, takes into account the enhancement to diffusion that arises from the effect of eddies increasing the length of tracer contours L and influencing the tracer gradients |∇c| through a combination of stirring and mixing. The magnitude of the effective diffusivity is thus directly related to the geometric complexity of tracer contours. Regions with relatively simple tracer contours and small κeff correspond to transport barriers, where the motion of Lagrangian particles is relatively regular. Generally, regions with a large κeff are representative of mixing regions where the motion of Lagrangian particles is typically chaotic.
c. Application to oceanic flows
In theory, the effective diffusivity for the surface ocean could be calculated from an observed tracer; however, the approach taken here is to calculate it from a tracer obtained by solving numerically the two-dimensional advection–diffusion equation (1), where u(λ, ϕ, t) is the nondivergent part of the near-surface geostrophic velocity field, as derived from altimetry. The infrastructure of the Massachusetts Institute of Technology (MIT) general circulation model (Marshall et al. 1997) was used to advect the tracer with an Adams–Bashforth time-stepping scheme and a simple centered second-order discretization in space that conserves c and c2 and introduces no spurious diffusion, as described in Marshall et al. (2006). Note that the effective diffusivity may also be calculated using a tracer field advected by velocity fields from numerical models (see, e.g., AMMS). We will present results for the Southern Ocean region; however, the method described could equally be applied to other oceans.3
A time series of sea surface height (SSH) anomalies was taken from the final combined processing of the Ocean Topography Experiment (TOPEX)/Poseidon (T/P) and European Remote Sensing Satellite-1/2 (ERS-1/2) data at 10-day intervals on a ¼° latitude × ¼° longitude grid, and the height relative to the geoid was obtained by subtracting the geoid model described in Lemoine et al. (1997). Figure 1a shows the SSH anomalies for 16 October 1998. There are distinct regions of strong variability in the SSH anomalies in the latitudes of the ACC and, in particular, in the regions of the East Australia Current, Southwest Indian Ridge, Brazil–Malvinas confluence, and the Agulhas Retroflection. These regions of large rms SSH anomalies correspond to regions of strong eddy kinetic energy (see Fig. 6 later).
To calculate the effective diffusivity we initialized a tracer field in the region of the Southern Ocean, poleward of 30°S, with the geographic distribution of the mean streamfunction, rescaled to have a value of 0 toward the north and 1 toward the south of the region [Fig. 2a shows a zoom of the tracer field at the start of an integration (30 December 1996) for a patch of the Pacific from 38° to 59°S, 160° to 140°W; the full Southern Ocean can be seen in Fig. 1 of Marshall et al. (2006)]. Transport in the cross-stream direction is typically slower than in the along-stream direction and, therefore, a tracer is expected in a first stage to adjust to broadly align with the mean streamlines and then subsequently to adjust more slowly across the mean streamlines. [Note that for a steady two-dimensional flow, such a two-stage tracer evolution can be described precisely (Rhines and Young 1983).] It can be seen from Figs. 2a–c that the tracer does, indeed, follow such an evolution. The first stage of the tracer evolution requires diffusive decay of the structure that arises purely from the initial condition. An initial condition aligned with the mean streamlines is chosen to minimize this adjustment time.
The tracer was evolved for one year4 with a value of numerical diffusion k chosen to prevent unrealistic structure on the grid scale and then output from the simulation at 10-day intervals. For each output the effective diffusivity κeff was calculated as a function of ϕe according to the method outlined in the appendix. The results are presented in Fig. 1b. Each tracer contour c = C corresponds to a particular equivalent latitude ϕe. Consequently, there is a value of κeff that corresponds to each tracer contour. To create the plot in Fig. 1b, we have colored each tracer contour according to its κeff values. It is important to realize that the κeff represents an average around a tracer contour, so care must be taken in interpreting the regional information in this plot [see Shuckburgh et al. (2009) for further discussion].
The instantaneous streamlines are overlaid in Fig. 1b. We have chosen to set the zero streamline to be approximately coincident with the ACC; positive streamlines lie poleward and negative equatorward following the conventions of Karsten and Marshall (2002). It can be seen that in the region of the ACC between streamline contours 0 and 6 × 104 m2 s−1 there are very low values of κeff (<1000 m2 s−1), equatorward of this up to the −5 × 104 m2 s−1 streamline are moderate values of κeff, and farther equatorward still higher values (>2000 m2 s−1). This pattern is quite different from that of estimates of the eddy diffusivity derived directly from altimetry data using the approach of Holloway (1986) (see also Keffer and Holloway 1988; Stammer 1998). Such calculations give eddy diffusivities that are typically enhanced in the regions of strong eddy kinetic energy.
As noted already, the assumptions behind the general approach of using large-scale velocities to advect tracers (or Lagrangian particles) have been discussed many times previously in the case of the atmosphere. The technique is, however, less well established for the ocean and is investigated by considering the effect of subsampling the advecting velocities to lower spatial resolution. For flows in which stratification and rotation are important, the small-scale tracer features are often governed predominantly by the large-scale flow (or in other words these flows often lie within the chaotic advection regime). Low-resolution velocity fields are expected to produce reliable estimates of small-scale tracer structures when there is a steep decay of the kinetic energy spectra, E(k) ∼ k−3 or steeper. As discussed by Waugh et al. (2006), there is some evidence that the energy spectra for the surface ocean fall within this regime [note, however, the recent results of Le Traon et al. (2008), which suggest a k−5/3 spectral slope]. In addition, Beron-Vera et al. (2008) find that they are able to reproduce the trajectory of observed drifters in the South Atlantic using the low-resolution velocities derived from altimetry. Returning to Fig. 2c, it can be seen that the tracer field exhibits smaller-scale structures than the streamfunction. The nominal resolution of the velocity field in this calculation is ¼°; however, the smaller scales may not be well represented in the satellite data. Eddy features can clearly be seen throughout the flow. If this streamfunction is subsampled to lower resolution by removing the highest wavenumbers, then many of the largest coherent structures are preserved at 1° resolution (Fig. 3a), but they are not preserved at 2° resolution (Fig. 3b). Correspondingly, there is considerable similarity between even the small-scale features of the tracer field after advection for one year with the full-resolution velocities (Fig. 2) and the 1° resolution velocities (Fig. 3a), but less similarity with the tracer field after advection with the 2° velocities (Fig. 3b). The equivalent length (Fig. 3c) is virtually identical even after one year of integration with velocities at ¼° and 1°, and there is even a strong similarity with velocities at 2°.
These results provide some confidence that it is, indeed, the larger-scale features in the flow that are responsible for creating the small-scale features in the tracer and, thus, that the flow falls within the chaotic advection paradigm. They also indicate that small-scale noise in the altimeter data is not having a strong influence on the results.
3. Robustness of the effective diffusivity diagnostic
The novel aspect of the effective diffusivity diagnostic is that the properties of the flow are inferred from a single advected tracer field with no direct reference to the flow field. There is, therefore, a concern that the diagnostic may be sensitive to details of the properties of the tracer. For steady flows, a tracer will (after the first stage of evolution) be aligned with the streamlines and hence the effective diffusivity can be written alternatively in terms of the streamfunction with no explicit reference to the tracer field (Rhines and Young 1983). For unsteady flows, Haynes and Shuckburgh (2000a,b) presented results from a number of simulations that indicated, in the atmospheric case, the effective diffusivity is determined primarily by the flow and is not strongly dependent on the properties of the particular tracer field. We now reexamine this issue for the ocean flow.
a. Dependence on initial conditions
We first examine whether the tracer indeed adjusts to a configuration for which the effective diffusivity is independent of the initial conditions, as argued above. Two tracer simulations were carried out: one was started 40 days after the other (30 December 1996 and 8 February 1997). Both simulations were at resolution in latitude and longitude with a numerical diffusivity of k = 50 m2 s−1. Figure 4a shows the effective diffusivity as a function of equivalent latitude for the tracer from the two simulations. A 50-day running mean in time and a ¼° running mean in equivalent latitude have been applied to the effective diffusivity results to remove some of the high-frequency noise arising from the calculations. High-frequency temporal variability is implausible and this filtering reduces the average difference between the effective diffusivity at one output time and the next to below 5%. The effective diffusivity distributions in the two simulations are shown to be very similar for each time plotted in Fig. 4a: 100, 190, and 280 days after the start of the second simulation (corresponding to 19 May, 17 August, and 15 November 1997). The latitudinal structure of the effective diffusivity is observed to be similar for each of these times, although there are some differences in the region from 50° to 40°S, which may indicate seasonal variability in this region; such seasonal variability will be explored further in section 5. The time evolution of the average (over all equivalent latitudes) percentage difference between the two simulations is shown in Fig. 4b. This confirms that after 100 days there is never more than about a 5% difference between the two simulations when the high-frequency filtering described above is applied to the results. On the basis of these results we argue that after a time of about 3 months the effective diffusivity is independent of the tracer initial condition. A similar result was found by Haynes and Shuckburgh (2000a), who argued for a 1-month “spinup” in the case of atmospheric flows. As noted by Kostrykin and Schmitz (2006) this spinup time will increase as the numerical diffusion k is decreased.
b. Dependence on resolution and numerical diffusion
The effective diffusivity is calculated from an instantaneous tracer field by evaluating Eq. (3). Some previous studies have considered artificial tracers that have been generated numerically by evolving a tracer field using a tracer transport model with velocities either based on observations (Haynes and Shuckburgh 2000a,b; Allen and Nakamura 2001; Shuckburgh et al. 2001; Marshall et al. 2006) or taken from a general circulation model (Nakamura and Ma 1997; Kostrykin and Schmitz 2006), and others have used tracer fields from observational data (Nakamura and Ma 1997). The terms Leq and Lmin in the effective diffusivity formula are straightforward to calculate from the tracer field; however, the diffusivity term k is more problematic. If k is unknown, but one assumes that it is a constant, then it is possible to determine the spatial structure of κeff and hence identify mixing regions and barriers in the flow. However, many studies require knowledge of the magnitude of the effective diffusivity and hence an estimate of k.
For an artificial tracer evolved using a numerical model, k is the numerical diffusion. This is either directly specified or its average value may be estimated by examining the decay with time of the globally averaged tracer variance (Allen and Nakamura 2001). Although Eq. (3) suggests at first sight a linear dependence of κeff on k, this is modulated by the influence of the diffusivity k on the growth of tracer contour length and hence on the term Leq. In the limit of large Péclet number, Pe = ULeddy/k, where U is the characteristic scale of the velocity, the diffusion is sufficiently small as to make its effects appreciable only on the small-scale tracer structures. We have demonstrated above that, for the ocean flows of interest, the small-scale tracer features are governed predominantly by the large-scale flow. Thus, for this class of flow, κeff may be expected to be determined primarily by the large-scale velocity field and to exhibit only weak dependence on k for large Péclet numbers [see Fig. 3 in Marshall et al. (2006) and discussions in Shuckburgh and Haynes (2003)]. The implication of this is that for the calculation of κeff, k can be chosen essentially arbitrarily as long as it is not so large as to dominate over the stirring effects of the large-scale advection nor so small as to prevent strong spatial variation of the tracer on the scale of the grid.
Figure 5 examines the effects of resolution and numerical diffusivity of the tracer integration on the effective diffusivity. The effective diffusivities for a series of patch calculations in the region 43°–53°S, 43°–53°E; [see the companion paper (Shuckburgh et al. 2009) for details of how these patch calculations are performed] are presented after one year of simulation with the filtering described above to remove high-frequency noise. The solid curves are the results from a tracer calculation with a resolution of in latitude and longitude, with numerical diffusivity of k = 1 m2 s−1 (black), k = 10 m2 s−1 (dark gray), and k = 100 m2 s−1 (light gray). In the weaker mixing region, between 48° and 42°S, there is little difference between the effective diffusivities arising from the k = 1 and 10 m2 s−1 calculations, whereas the effective diffusivity from the k = 100 m2 s−1 calculation is noticeably higher. This suggests that, for weaker mixing regions, the numerical diffusivity is dominating over the effects of advection at k = 100 m2 s−1. However, in the stronger mixing region around 50°S, the effective diffusivity from the k = 1 m2 s−1 calculation is noticeably smaller than the effective diffusivity arising from the calculations with higher numerical diffusion. This suggests that for the stronger mixing region, k = 1 m2 s−1 is not sufficient to prevent tracer structure arising on the grid scale, which is providing an additional diffusivity not accounted for here. Therefore, for these calculations at resolution, k = 10 m2 s−1 seems a good compromise.
The dashed curves in Fig. 5 are the results of a calculation with a resolution of in latitude and longitude, with numerical diffusivity as for the calculation. In the weak mixing region the results of the calculation are similar to those from the calculation. However, in the stronger mixing region between 50° and 52°S the results of the k = 1 and k = 10 m2 s−1 at resolution calculation are reduced compared to those of the resolution calculation. This suggests that a numerical diffusion of k = 10 m2 s−1 is not sufficient to present structure arising on the grid scale in this region at .
The above results are consistent with the conclusion of Marshall et al. (2006) that at resolution a numerical diffusivity of k = 50 m2 s−1 provided an effective diffusivity that was relatively independent of the numerical diffusivity: k = 100 m2 s−1 is too large in weak mixing regions and k = 10 m2 s−1 is too small in strong mixing regions.
c. Dependence on tracer properties
The results of this section have demonstrated that the effective diffusivity is not strongly dependent on the details of the properties of the tracer used in its calculation. We conclude that, as for atmospheric flows, the effective diffusivity may be considered a robust quantifier of the eddy diffusivity for oceanic flows.
4. Spatial distribution
The spatial structure of effective diffusivity, as presented in Fig. 1, has low values at the core of the ACC, with moderate values on its flanks, and the highest values in the most equatorward latitude bands (recall that the effective diffusivity is a streamwise-average quantity). In this section we attempt to explain this basic spatial structure by exploring how it relates to other measures of the flow. Figure 6 presents a number of such measures. Figure 6a provides an estimate of the phase speed (c) of the eddies, calculated by C. Hughes from observational data using Radon transforms [the same data have been presented previously by Smith and Marshall (2009)]. Figures 6b and 6c indicate the time-mean eddy kinetic energy (EKE) and time-mean zonal current (u) from altimetric observations, averaged over the time period from January 1996 to December 2000 (the data have been smoothed with a 2° running mean).
Throughout much of the Southern Ocean, away from the ACC where the zonal currents are weak, the sea surface height shows westward propagating anomalies with phase speeds less than −2 cm s−1. In the ACC the anomalies propagate eastward, but the phase speeds (2–5 cm s−1) are less than the zonal current, which can exceed 10 cm s−1. Linear baroclinic instability theory (Green 1970; Marshall 1981; Killworth 1997) suggests that the eddy diffusivity should scale with the eddy kinetic energy, zonal current, and eddy phase speed according to EKE/(u − c)2. This leads to an expectation of mixing and enhanced eddy diffusivities in the vicinity of critical layers where u = c and reduced diffusivities in regions such as the ACC where the zonal flow significantly exceeds the phase speed of the eddies (Smith and Marshall 2009). This picture is consistent with the potential vorticity (PV) distribution. As noted by Marshall et al. (2006), weak PV gradients are observed in the mixing region equatorward of the ACC and strong gradients (typically associated with eddy transport barriers) are found roughly coincident with the core of the ACC.
Figure 6d presents the spatial distribution of EKE/(u − c)2, averaged over the period January 1996–December 2000, calculated from the smoothed data in Figs. 6a–c. Significant longitudinal variability exists. Broadly, however, the latitudes of the ACC are characterized by lower values of this quantity than the equatorward latitudes, in agreement with the meridional variation in effective diffusivity (see Fig. 1). This latitudinal dependency can be seen more clearly in Fig. 7, which compares the effective diffusivity as a function of equivalent latitude (black) with streamwise-average EKE/(u − c)2 (gray). The structure of lower values in the latitudes of the ACC, and higher values equatorward of this, is clearly evident in both quantities.
Further evidence for this dependency of the eddy diffusivity on the flow parameters is found in the calculations of AMMS. In calculations based on model velocities constrained by observations, they find enhanced eddy diffusivities below the surface in the Southern Ocean, at the steering level of baroclinic waves. The independent studies of Smith and Marshall (2009) also find a similar dependency. Thus, from a number of different studies, a collective picture is emerging: consistent with the theory of baroclinic instability and potential vorticity mixing, of reduced eddy diffusivity at the core of the ACC, and enhanced diffusivity at critical layers on its flanks.
5. Temporal variability
a. Temporal variability of effective diffusivity
To investigate the temporal variability of eddy diffusivity, we consider two dates six months apart: 16 October 1998 (Fig. 1) and 9 April 1999 (Fig. 8). Comparison of the SSH anomalies for the two dates indicates (i) the sea surface is elevated over much of the Southern Ocean in April 1999 (see comments below) and (ii) the eddy activity is slightly enhanced in April 1999 compared with October 1998. Consistent with this enhanced eddy activity, the eddy diffusivity in the most equatorward band is seen to be enhanced in April 1999. Interestingly, however, the eddy diffusivity in the ACC band is reduced in April 1999 compared with October 1998. Analysis of the flow indicates that the mean zonal current in the ACC band was weak in October 1998 (7.8 cm s−1, just over one standard deviation below the time mean for the period Janurary 1996–December 2000), whereas it was very strong in April 1999 (8.3 cm s−1, nearly two standard deviations above the time mean). It was noted above that the eddy diffusivity scales with EKE/(u − c)2, so the increased ACC zonal current in April 1999 would be expected to reduce the eddy diffusivity, as is observed. It should be noted that the increase in eddy activity and the increase in speed of the zonal current in April 1999 may be dynamically linked through wave–mean flow interactions, and the entire pattern of variability in flow parameters and eddy diffusivity is consistent with the expectations based on potential vorticity mixing (Marshall et al. 2006; Haynes et al. 2007; Dritschel and McIntyre 2008).
The temporal variability of the effective diffusivity is further investigated in Fig. 9. The κeff results from a series of nine 12-month tracer simulations are presented. Each of the dark gray curves represents a simulation initialized on 30 December; the light gray curves on 8 July. The lowermost curve, with an average value of about 850 m2 s−1, is an average between streamlines in the region of the ACC [with values 0 and 6 (×104 m2 s−1), as marked in Fig. 1, corresponding to equivalent latitudes 49°–56°S]. The middle curve with an average value of ∼1450 m2 s−1 corresponds to the region just equatorward of this region on the flanks of the ACC (with streamline values from 0 to −5 × 104 m2 s−1, equivalent latitudes 41°–49°S). The top curve with an average value of 2300 m2 s−1 corresponds to streamlines from −5 to −8 (×104 m2 s−1) and equivalent latitudes 33°–41°S.
The 3-month spinup period for each calculation is not plotted, and the overlap period between the dark and light gray curves is plotted as thin contours. For the most part, the agreement between the two calculations during the overlap period is very good, confirming the robustness of the calculation. The exception is the equatorwardmost region for the calculations initialized on 30 December 1997, 30 December 1998, 8 July 1999, and 30 December 1999. It is thought that the longer spinup required at these times in this region is because the variability in this region means that it takes longer for the tracer to adjust from the initial conditions of the time-mean streamfunction during the first stage of evolution: initializing instead with the instantaneous streamlines might help to reduce this adjustment time.
The black curve shows the results after they have been subjected to a fifth-order Butterworth low-pass filter with a 90-day cutoff frequency to focus on the seasonal variations. There is no obvious trend in any of the regions. The greatest variability is observed in the equatorwardmost band, although the variability in the three bands is approximately the same if they are normalized by their average value. In the region of the ACC (lowermost curve), the strongest barrier is observed from January to July 1998 and from January to July 2000, while the weakest barrier is observed from August 1998 to February 1999 and August to November 1999. In the equatorwardmost band (highest curve), enhanced mixing is predominantly in fall: the greatest mixing is observed from January to August 1999, with smaller peaks in August–September 1998, March–May 1998, May–July 1997, and November 1999–May 2000.
The seasonal cycle of eddy diffusivity is further investigated in Fig. 10. Equatorward of the ACC (Fig. 10, top panel) the eddy diffusivity appears to be greatest in fall (March–April) and spring (November) and least in winter (June–October). The eddy kinetic energy follows a similar seasonal cycle in this region. There is some anticorrelation with the eddy diffusivity in the ACC region (Fig. 10, bottom panel), being least in fall (February–June) and greatest in winter/spring (July–December). This weak anticorrelation might be a consequence of the stronger mixing equatorward of the ACC acting to sharpen the jets of the ACC as a consequence of its action on potential vorticity gradients, thus creating a stronger barrier at the core of the ACC (Dritschel and McIntyre 2008; McIntyre 2008).
b. Temporal variability of eddy activity
To investigate the causes of the temporal variability in eddy diffusivity we look first in detail at the time evolution of the eddy activity by considering the eddy kinetic energy averaged over the equivalent latitude bands used in Fig. 9, EKE. We also consider a quantity based on the instantaneous sea surface height anomalies
where h′ is the SSH anomaly and α is a constant of proportionality, which we set to 0.26 (with the average defined as before). This quantity is related to the eddy diffusivity introduced by Holloway (1986), which is based on the rms h′. Figures 11 and 12 present the time evolution of EKE and κalt, respectively, over the time period 1996–2001. There are three clear peaks in κalt in fall 1998, 1999, and 2000 in the bands equatorward of the ACC; however, the only significant peak in EKE in these regions is in fall 1999. The seasonal heat storage in the ocean, yielding thermosteric height changes, appears to be the leading contribution to SSH variance5 (and hence to κalt) equatorward of 50°S; poleward of this and in the deep basins, the barotropic response to wind forcing dominates the large-scale SSH variability (Vivier et al. 2005). This large-scale variability likely dominates the temporal variability in Fig. 12, but probably has rather little impact on the EKE (or eddy diffusivity). This means that care needs to be taken in interpreting diffusivity measures based on SSH anomalies arising directly from the altimeter data without taking into account such large-scale variability.
Meredith and Hogg (2006) investigated the temporal variability of annual-mean EKE (derived from altimeter measurements) averaged over boxes in the Pacific, Atlantic, and Indian Ocean sectors of the Southern Ocean. The boxes were constructed to enclose the high EKE areas of the northern ACC, but exclude regions where the ACC interacts with noncircumpolar current systems such as the Brazil and Agulhas Currents. They considered the time period 1993–2004 and noted a strong increase in annual-mean eddy kinetic energy in the Pacific and Indian sectors in 2000, 2001, and 2002 and in the Atlantic sector in 2001. Examination of the latitude–longitude distributions of SSH anomalies and EKE reveals that, within the ACC region, the greatest changes are located in a few specific regions, such as the Eltanin and Udinstev Fracture Zones and the Campbell Plateau in the Pacific (see Fig. 2 in Meredith and Hogg 2006). In Fig. 11 we have averaged the EKE around all longitudes over three equivalent latitude bands. Given the regional variation of EKE, this leads to some differences compared to the results of Meredith and Hogg. The most noticeable is that the annual-mean EKE (dashed line) is found to be as high in 1999 as in 2000 (with both years being higher than 1996–98). Inspection of the latitude–longitude distributions indicates a significant difference in EKE in the Pacific Fracture Zone between the two years, which contributes greatly to the Pacific box of Meredith and Hogg (2006), but less so to the larger averaging region in Fig. 11.
c. Temporal variability of flow parameters
We noted in section 4 that linear baroclinic instability theory suggests that the eddy diffusivity should scale with EKE/(u − c)2 and found the spatial distribution of these two quantities to be in relatively good agreement. To investigate this relationship further, we consider in Fig. 13 the temporal variations of the two quantities. (For clarity of presentation, we have applied a fifth-order Butterworth low-pass filter with a 180-day cutoff frequency; however, the results are not strongly dependent on the choice of cutoff frequency.)
In the region of the ACC (Fig. 13, bottom curve), there is a statistically significant correlation between κeff (black) and EKE/(u − c)2 (dark gray) with a value of r = 0.56. In this region variations in the mean flow dominate variations in the eddy diffusivity, with a statistically significant correlation between κeff and u of r = −0.47 (for κeff and EKE the correlation is also statistically significant in this region, but r = 0.30). In contrast, in the equatorwardmost band variations in the eddy kinetic energy dominate, with a statistically significant correlation between κeff and EKE of r = 0.69 (the values of u and c in this band may be complicated by the inclusion of some basin-scale barotropic modes). Thus, the mean flow and interactions between the eddies and the mean flow play important roles in determining the temporal as well as the spatial variability in eddy diffusivity.
In their investigation of the temporal variability of EKE, Meredith and Hogg (2006) suggested that the increase that they observed in EKE in 2000–02 may have been a 2–3-yr lagged response to a peak in the circumpolar eastward wind stress associated with the Southern Annular Mode [SAM, the leading mode of variability in the Southern Hemisphere, characterized by geopotential height perturbations throughout the troposphere and stratosphere of opposing signs in the polar cap region and the midlatitudes; Thompson and Wallace (2000)]. We have included the 2-yr lagged SAM index, obtained by taking the emperical orthogonal functions of monthly mean 1000-hPa height anomalies poleward of 20° latitude, in Fig. 13 (light gray dash). With the rather short time period under consideration, we do not find a significant correlation between the 2-yr lagged SAM and the effective diffusivity in any of the latitude bands (we also investigated other lags from 0 to 3 yr but similarly found no significant correlation).
6. Conclusions and discussion
The purpose of this paper was to verify the robustness of the effective diffusivity diagnostic in an oceanographic context. We have chosen to use the Southern Ocean, a region of strong eddy activity, as the focus for our analysis.
Detailed analyses have shown that the effective diffusivity diagnostic, as implemented by Marshall et al. (2006), is indeed a robust diagnostic of the eddy diffusivity. It has been shown that the surface velocities fall into the chaotic advection regime and, therefore, small-scale inaccuracies in the satellite altimetry do not greatly impact the effective diffusivity calculation. (Note that errors in the mean flow may have a greater impact; this will be investigated in a future study.) It has also been shown that the effective diffusivity is not sensitive to the initial conditions of the tracer used in the calculation after a spinup time of approximately 3 months. As discussed by Marshall et al. (2006), the numerical diffusivity used in the calculation should be chosen such that it is small enough for the effective diffusivity to be largely independent of it, but large enough to prevent spurious small-scale tracer structure in mixing regions.
The differences between eddy features in the atmosphere and in the ocean were described in section 2. In the atmosphere strong barriers to isentropic eddy transport exist at the edge of the stratospheric polar vortex and at the subtropical jet in the tropopause region. Calculations of equivalent length give values of close to Lmin for the barrier at the polar vortex edge (Haynes and Shuckburgh 2000a) and approximately 1.5 × Lmin for the subtropical jet barrier (Haynes and Shuckburgh 2000b) in winter. It can be seen in Fig. 3 that the minimum value of equivalent length reached in the ACC region is approximately 4 × Lmin. It is therefore perhaps more appropriate to consider the ACC to be a region of weaker mixing, rather than a barrier to eddy transport.
It has been shown that the spatial and temporal variability in eddy diffusivity follows that anticipated from linear baroclinic instability theory, broadly following the streamwise-average EKE/(u − c)2. This accounts for the low values of eddy diffusivity in the vicinity of the ACC. Although the eddy kinetic energy is high throughout much of the ACC region, the phase speed of the eddies is significantly slower than the zonal flow, resulting in a reduced eddy diffusivity. The consequence of this relationship is that both the eddy kinetic energy and the zonal current need to be taken into account in characterizing the eddy diffusivity.
Analysis of the temporal variability of eddy diffusivity indicates significant variability on seasonal time scales with a magnitude of ±10%–15% of the mean. Equatorward of the ACC there appears to be a tendency for enhanced values of eddy diffusivity in fall, but there is no indication of a systematic trend over the five years of our calculation. Over this time period, 1996–2001, observations (Meredith and Hogg 2006) indicate a substantial peak in annual-mean eddy kinetic energy at the end of the period (2000), at least in the Pacific and Indian Ocean sections. Our calculations suggest that the response in terms of the eddy diffusivity may be less remarkable. The residual-mean framework of Marshall and Radko (2003) relates the eddy diffusivity to the strength of the overturning circulation and to the baroclinic transport of the ACC. The lack of a clear direct relationship between the eddy kinetic energy and the eddy diffusivity suggests that one might not expect a direct relationship between variations in the eddy kinetic energy and variations in the strength of these dynamical quantities. This has implications for the effect on the dynamics of the Southern Ocean of changes in the wind stress (see also Hallberg and Gnandesikan 2001; Hogg and Blundell 2006; Hogg et al. 2008) and the role of eddies in the response. This is a topic of much debate since observations indicate the wind stress is increasing (Thompson and Solomon 2002) and there are potential implications for the uptake of carbon dioxide in the Southern Ocean (Le Quéré et al. 2007). Hallberg and Gnandesikan (2006) find that, in models with explicitly resolved eddies, the response of the Southern Ocean meridional overturning circulation to changes in the wind stresses is significantly smaller than in models with parameterized eddies, as the eddy-driven circulation compensates for changes in the Ekman transport. In addition, the observational record shows only small changes (5% on interannual time scales) in ACC transport despite changing winds; model results (e.g., Hallberg and Gnandesikan 2006) also indicate a weak ACC transport response. One of the advantages of the effective diffusivity diagnostic is that it can be applied to both models and observations. A future study will directly investigate the influence of changes to the wind stress on the eddy diffusivity and the eddy-driven overturning circulation in a model context.
EFS would like to thank Mike Meredith for helpful discussions. EFS is supported by a Natural Environment Research Council fellowship. JM acknowledges the support of the National Science Foundation (Polar Program). Two anonymous reviewers made very useful comments on the manuscript.
Effective Diffusivity Calculation
The effective diffusivity, κeff, is calculated as a function of the equivalent latitude ϕe from the simulated tracer field as follows.
(i) The correspondence between tracer contours and equivalent latitudes is tabulated—landmasses in the region are treated by assigning them values below (above) the minimum (maximum) tracer and thus removing them from the domain of calculation. The land of Africa, South America, and Australia–New Zealand is assigned a value less than 0 and that of Antarctica a value greater than 1.
(ii) The equivalent length Leq is calculated as defined in Eq. (4)—the gradient of the tracer is calculated at each grid point, taking care at the land–sea boundaries, and its square is integrated over the area bounded by the desired tracer contour and then differentiated with respect to area by taking finite differences. The resulting quantity is then divided by the square of the areal gradient of the tracer.
(iii) The minimum length Lmin is calculated—this is estimated numerically for the region of interest by calculating Leq for a tracer that has been evolved to a steady state in terms of contour lengths under the influence of excessive diffusion (e.g., k = 104 m2 s−1).
(iv) The total numerical diffusivity is estimated, which will include both the explicit numerical diffusivity k and a contribution (which may be negligible) resulting from the numerical advection scheme—this is achieved by examining the decay with time of the globally averaged tracer variance.
(v) The effective diffusivity is constructed by multiplying the total numerical diffusivity by the equivalent length and normalizing by the minimum length.
Corresponding author address: Emily Shuckburgh, British Antarctic Survey, High Cross, Madingley Road, Cambridge, CB3 0ET, United Kingdom. Email: firstname.lastname@example.org
Here Leddy is smaller than the Rhines scale, approximately 100 km in this region, at which one expects the energy cascade in geostrophic turbulence from smaller to larger scales to be arrested.
Scaling arguments can be used (Vallis 2006) to relate the magnitude of the eddy kinetic energy to the magnitude of the mean kinetic energy multiplied by (Leddy/LD)2.
Note, however, that the assumptions of geostrophy break down at low latitudes and the estimates of sea surface variability are invalid over sea ice at high latitudes.
Simulations were halted after one year since the numerical accuracy of the calculation of effective diffusivity relies on a strong latitudinal gradient of tracer.