Hughes et al. (2010)

Title:

Identification of jets and mixing barriers from sea level and vorticity measurements using simple statistics

Key Points:

  • Skewness and kurtosis in sea level data can be used to identify oceanic fronts and mixing barriers.

  • Regions indicative of strong jets are associated with zero skewness and low kurtosis.

  • Geostrophic relative vorticity was the preferred variable to represent local dynamics.

Keywords:

skewness, kurtosis, sea level, vorticity, jet, meander, mixing, AVISO, altimetry

Corresponding author:

Chris W. Hughes

Citation:

Hughes, C. W., Thompson, A. F., & Wilson, C. (2010). Identification of jets and mixing barriers from sea level and vorticity measurements using simple statistics. Ocean Modelling, 32(1–2), 44–57. doi:10.1016/j.ocemod.2009.10.004

URL:

https://doi.org/10.1016/j.ocemod.2009.10.004

Abstract

The probability density functions (PDFs) of sea level and geostrophic relative vorticity are examined using satellite altimeter data. It is shown that departures from a Gaussian distribution can generally be repre.sented by two functions, and that the spatial distribution of these two functions is closely linked to the skewness and kurtosis of the PDF. The patterns indicate that strong jets tend to be identified by a zero contour in skewness coinciding with a low value of kurtosis. A simple model of the statistics of a mean.dering frontal region is presented which reproduces these features. Comparisons with mean currents and sea surface temperature gradients confirm the identification of these features as jets, and confirm the existence of several Southern Ocean jets unresolved by drifter data. Diagnostics from a range of idealized eddying model simulations show that there is a strong, simple relationship between kurtosis of potential vorticity and effective diffusivity. This suggests that kurtosis may provide a simple method of mapping mixing barriers in the ocean.

Introduction

More than a decade of precise satellite altimeter measurements has produced long time series of sea level h, over most of the global ocean. As a way of representing the different variations seen in different regions, it is common to plot the root-mean-square devi.ation of sea level from its mean, r.h., or the geostrophic eddy kinetic energy (EKE), large values of which tend to be associated with strong, eddying currents. More recently, Thompson and Demirov (2006) noted that there is further interesting information in the skewness $S(h)$ of the sea level distribution, which tends to be positive on the poleward side of strong eastward currents, and negative on the equatorward side. They interpreted this as indica.tive of the intermittent passage of warm core eddies or warm meanders across points poleward of a jet, and of cold core eddies or meanders on the equatorward side.

This prompts the question of whether higher moments of variability contain further information, and how those moments might be interpreted. More generally, what is the spatial variation in the PDF of sea level, and what causes that variation? In considering this question, we found that sea level, being a non-local variable from a dynamical point of view, is sometimes more difficult to interpret than the more localized relative vorticity, so in much of the work presented here we focus on the PDF of relative vorticity $zeta$.

We show below that, in fact, maps of skewness $S$ and kurtosis $K$ contain most of the information which can be reliably extracted from PDFs at each position. These are defined as $$ S(chi) = frac{frac{1}{n}sum_nchi^3}{(frac{1}{n}sum_nchi^2)^{frac{3}{2}}} $$, (1) $$ K(chi) = frac{frac{1}{n}sum_nchi^4}{(frac{1}{n}sum_nchi^2)^2} $$, (2) P3=2 ; 1 x2 n Pn

where $chi$ is the deviation of some quantity from its mean value. Skewness is a measure of asymmetry of the PDF. Where rare, large events predominantly have one sign, this is the sign of the skewness; a Gaussian distribution is symmetric and therefore has a skewness of zero. Kurtosis is traditionally described as a measure of ‘‘peakiness” of the PDF. High kurtosis means that rare, large amplitude events occur more frequently than would be expected for a Gaussian distribution, producing a more sharply peaked PDF with broader tails than a Gaussian. Low kurtosis means that events much larger than typical are rarer than would be expected from a Gaussian distribution, and the associated PDF tends to be flat near its peak, or even bimodal, and drops to small values faster than a Gaussian. The smallest possible kurtosis is 1, and this results from a variable taking only two possible values, with equal probability. A Gaussian distribution leads to a kurtosis of 3. For this reason, the quantity $K - 3$ is often used, sometimes termed excess ‘‘kurtosis” and sometimes simply ‘‘kurtosis” . We will use the definition given above and, when 3 is subtracted, this will be stated explicitly.

Our dataset consists of the AVISO gridded sea level anomaly (corrected for tides and inverse barometer effect) supplied weekly on a 1/3 degree Mercator grid. While the grid is 1/3 degree, the spacing between altimeter tracks is 0.72 degrees of longitude for ERS/Envisat, and 2.78 degrees for Topex or Jason, so we should not expect to resolve all features less than about 100 km across. We use the ‘‘reference” version of the data in which sampling is always from two altimeters (Topex/Poseidon or Jason 1 and ERS or Envisat) in the same two orbits, and consider 690 weeks spanning the period 5th April 1995–11th June 2008 which follows the period when ERS-1 was not in the standard 35-day repeat orbit. Geostrophic relative vorticity anomalies are calculated by first-order centered differencing of the sea level to give geostrophic velocity, and the same form of differencing to then give vorticity at the same grid points as the sea level data. These are divided by the local Coriolis parameter $f$ so that positive values represent cyclonic, and negative values anti-cyclonic, relative vorticity.

As we are interested in the shape of the PDFs rather than their widths, we consider only time series which have had their means removed and have been divided by their standard deviations. In fact, we also subtract an annual and semiannual cycle and a linear trend, simultaneously Þt to each time series. It is also worth noting a limitation of the geostrophic relative vorticity calculation. The effect of nonlinear terms in the momentum equation is such that the geostrophic velocity (and hence vorticity) is an underestimate in cyclonic vortices, and an overestimate in anticyclonic vortices. This potential bias should be kept in mind, although it is not particularly important for the main issue we consider here, which focuses on the step-like structure in jets.

Binning the resulting normalized PDFs from each grid point into bins 0.01 standard deviations wide, and averaging over the globe (excluding the band within 5 degrees of the equator in the case of vorticity, to avoid the equatorial singularity), produces the average PDFs shown on both logarithmic and linear scales in Fig. 1. The dashed lines show the corresponding Gaussian PDF for reference. We see that the average shape of the PDF for both sea level and relative vorticity is fairly close to Gaussian, but significantly above Gaussian in the wings and centre of the distribution (and below Gaussian in between). Gille and Sura (submitted for publication) look at the overall global relationship between skewness and kurtosis seen from such PDFs, and similar characterization of the velocity PDFs from altimetry has been investigated by Schorghofer and Gille (2002) and Gille and Llewellyn-Smith (2000). This average shape, however, conceals a large geographic variation. Characterization of this spatial variability is the subject of the next section.

Fig. 1: Global average, normalized PDF of (left) sea level, and (right) geostrophic relative vorticity $zeta/f$f from gridded altimetry. The same data are shown on (top) logarithmic and (bottom) linear scales. In each panel, the dotted line is the corresponding Gaussian PDF.

Spatial variations of the PDFs

Skewness and kurtosis are the first two in an infinite series of moments which characterize the shape of the PDF, following the mean and standard deviation which characterize its position and width. However, it is not immediately clear that these are the parameters best suited to describing the spatial variations in shape. Furthermore, as kurtosis is particularly sensitive to rogue large values in a dataset (it is sometimes used as an indicator of bad data), it is quite possible for much of the shape of the PDF to be invisible to the kurtosis in cases where one particularly large value occurs. For this reason, we seek to characterize the typical variations in shape of the PDF, ignoring the very extreme wings of the distribution.

To do this, we use empirical orthogonal function (EOF) analysis. This is often used to calculate spatial patterns associated with the principal modes of variability of time series given at each grid point. Here, instead of time, we are using the x-axis of the normalized PDF graph at each grid point. Rather than calculate EOFs of the complete PDF, we choose fist to subtract the Gaussian PDF so that we will clearly be mapping departures from Gaussian behavior. We obtain a set of orthonormal eigenvectors, the first of which explains the largest possible amount of spatial variability in the PDF, the second of which describes as much as possible of the remaining variability after removing the effect of the first mode, and so on. The associated eigenvalues give the amount of variance explained by each mode.

There are not enough points in a time series at one point to produce a PDF with such high resolution as in Fig. 1, so here we increase the bin size to 0.1 standard deviations, and consider only the 80 points between -4 and +4 standard deviations (this avoids the inclusion of any very extreme values). The analysis was also repeated with a bin size of 0.4 standard deviations, and only 20 points per PDF, with very similar results to those presented here.

The first 15 eigenvalues of the resulting EOFs are plotted in Fig. 2, as percentages of the total variance (diamonds for the sea level EOFs and crosses for relative vorticity). In both cases, the first two eigenvalues stand clearly above the background continuum, together explaining 28.8% of the sea level variance and 22.7% of relative vorticity variance. In the case of sea level, the second two EOFs may also plausibly be said to stand above the background. Inspection of the associated eigenvectors clearly shows the remaining EOFs to be associated with statistical noise, representing simply the exchange of probability between neighbouring bins (using coarser bins reduces the amplitude of this noise, leaving the first two EOFs to account for 61.7% of variance for sea level and 60.6% for relative vorticity, but does not lead to identification of any more significant modes).

Fig. 2: Eigenvalues of the first 15 EOFs of the PDFs for (diamonds) sea level and (crosses) relative vorticity. Eigenvalues are normalized so that they sum to 100, and each represents the percentage of variance explained by the corresponding eigenvector.

Since we are interested in relating the variations in PDFs to skewness and kurtosis, it would be nice if the EOFs fell naturally into symmetric and antisymmetric categories. They do not, but the importance of symmetry is clear from the fact that very nearly symmetric and antisymmetric functions can be formed from the obvious pairs of EOFs. If we take the first two EOFs of sea level, we can construct from these by a rotation, two different orthonormal vectors $mathbf{C}$ and $mathbf{D}$ chosen to minimize the quantity $(mathbf{C} +mathbf{C}^r)^2 + (mathbf{D} -mathbf{D}^r)^2$, where the superscript $r$ represents the reverse of the vector (i.e. the components listed in reverse order). This minimizes the sum of the squares of the symmetrical component of $mathbf{C}$ and the antisymmetric component of $mathbf{D}$. In other words, $mathbf{C}$ represents a linear combination of the first two EOFs which is close to antisymmetric, and $mathbf{D}$ an orthogonal linear combination which is close to symmetric. The results of such a minimization applied to the pairs of EOFs identified from Fig. 2 are shown in Fig. 3. The solid lines are constructed from pairs of EOFs of the sea level PDFs, and dashed lines from relative vorticity. To show how closely these approximate the ideal symmetry, crosses represent the exactly antisymmetric or exactly symmetric components of the corresponding sea level EOF curves. Success in this minimization is not inevitable. For example, no combination of the first and third EOFs can produce anything like a symmetric and antisymmetric pair of functions. This success suggests that it is meaningful to combine the EOFs in this way.

Fig. 3: Rotated eigenvectors of the first 4 EOFs (solid) of the PDFs for sea level and the first 2 EOFs (dashed) for relative vorticity. Crosses are half of the corresponding sea level vector plus or minus its reverse, illustrating the corresponding perfectly symmetrical or antisymmetric curves.

The projection of the PDFs onto these rotated EOFs is shown in Figs. 4 and 5. The pattern of the first rotated EOF for sea level is clearly very similar to the skewness map presented by Thompson and Demirov (2006). Major eastward currents are delineated by the boundary between large positive and large negative values, and preferred paths of eddies, such as in the Alaskan Stream along the Aleutian island chain at the northern boundary of the PaciÞc (Ueno et al., 2009), and southwest of Australia (Morrow et al., 2004) also stand out. The tropical PaciÞc is also prominent. As Thompson and Demirov (2006) point out, this is the result of the dominance of the large 1997Ð1998 El Nino event in the region. This highlights a disadvantage of considering sea level: large scale sea level signals can dominate the time series at a given point, but there may be little or no associated current or other dynamical change at that point. Dynamically, sea level is a non-local variable on time scales of several days or longer. For this reason, the localization of physical processes is better determined using vorticity or potential vorticity. This is illustrated for relative vorticity in Fig. 5. The spatial pattern for the first rotated EOF is clearly related to that for sea level (with a sign change since a positive localized sea level anomaly represents an anticyclonic circulation), but the relative vorticity map shows finer scales and more localized structures (particularly in the Southern Ocean) than the corresponding map based on sea level. Using relative vorticity also avoids a dominant influence of the El Nino event, making it possible to pick out clear structures to the west of Hawaii and Central America.

Fig. 4: Projection of the sea level PDFs onto the first and second rotated eigenvectors from Fig. 3.

Fig. 5: Projection of the relative vorticity PDFs onto the first and second rotated eigenvectors from Fig. 3.

The projections of sea level and relative vorticity PDFs onto the second rotated EOFs bear a similar relationship to each other (without the sign change, since the second rotated EOF is close to symmetrical). Strong currents now stand out as low values of the projection, and these features are more clearly deÞned in the relative vorticity plot.

From the forms of the EOFs, we expect the projections to be closely related to skewness and kurtosis respectively, and Fig. 6 shows that this is the case for relative vorticity (the equivalent maps for sea level skewness and kurtosis are not shown, but bear a similar relation to the sea level EOF projections). The EOF projections can be described as slightly “cleaner”: definition of some features is better, such as the Loop Current in the Gulf of Mexico, and features along the Falklands shelf. There are also some isolated spots of large skewness and kurtosis which are not visible in the EOF maps, which suggests that they may be due to rogue data points (a few time series which have been checked are consistent with this, consisting of a rather ßat curve with one sudden spike). A problem with the gridding in the region north of Iceland and Norway also shows up much more clearly in skewness and kurtosis than in the EOF maps. However, it can generally be said that maps of skewness and kurtosis are capturing the same spatial variability as the two EOF maps and, since only two EOFs are above the background noise for relative vorticity, that means that skewness and kurtosis capture all that can currently be captured about spatial variations in the PDF of relative vorticity. In the case of sea level, the spatial patterns associated with the third and fourth rotated EOFs (not shown) appear very noisy except in the vicinity of the very strongest currents, suggesting that they do contain physically meaningful information, but only in these select regions.

Fig. 6: Skewness and kurtosis of relative vorticity time series.

It is difficult to assess the statistical noise levels on these estimates, as the weekly sea level maps are not independent. A standard result estimates the standard error on skewness as $sqrt{6/n}$, and that on kurtosis as $sqrt{24/n}$, where $n$ is the number of independent samples, although kurtosis has a strongly asymmetric distribution, having a minimum value of 1, Gaussian value of 3 and inÞnite maximum possible value for a continuous PDF. As a guide, the 5 and 95 percentiles for skewness and kurtosis generated from a Gaussian distribution sampled at a Þnite number $n$ of independent times have been calculated by Monte Carlo simulation (using 100,000 random input time series) for $n = 690$ (the number of points in the time series), $n = 300$, and $n = 100$. The results are shown in Table 1. The strongest justification for the significance of these results, however, is the spatial coherence and physical plausibility of the resulting patterns.

Table 1: Values of skewness ($S$), kurtosis minus 3 $(K - 3)$ and projections of rotated EOF1 and rotated EOF2 generated by chance from Þnite samples of length $n$ taken from Gaussian distributions. The 5th and 95th percentiles are given.

$n$

%

$S$

$K - 3$

EOF1

EOF2

690

5

-0.14

-0.36

-0.10

-0.11

95

0.14

0.16

0.10

0.07

300

5

-0.21

-0.48

-0.15

-0.16

95

0.21

0.29

0.15

0.12

100

5

-0.36

-0.73

-0.26

-0.27

95

0.36

0.51

0.26

0.23

Interpretation of low kurtosis regions

As stated in the introduction, the lowest possible kurtosis is 1, and occurs when the variable ($zeta/f$ or $h$) is always at one of two values, and occupies those two values with equal probability. Consider a sharp jet, which may be approximated as a step in sea level. If that jet meanders, then sea level within the region over which it meanders will always be at one of two values (either to the left, or to the right of the jet). Along some line, the jet will spend equal amounts of time to the left and to the right of that line. On that line, the kurtosis of sea level would be 1.

Off the line, the fraction of time spent in each state varies, becoming more asymmetrical with distance from the centre line until beyond the limit of meanders where sea level becomes constant. If the fractions of time spent at each value are $A$ at the lower value, and $B = 1 - A$ at the higher value, then it is simple to calculate that $K = -3 + 1/P$, $S = (A -B)/sqrt{P}$, $S^2 = -4 + 1/P$, where $P = AB$ has a maximum possible value of 1/4. It will be noticed that this rather extreme situation gives the relationship $K = S^2 + 1$. This in fact represents the minimum possible $K$ for a given $S$ (Pearson, 1916). Idealized though this model is, it does capture the main features surrounding jets in Fig. 6: low kurtosis and zero skewness at the centre (mean position) of the jet, with kurtosis growing to large positive values on either side, as skewness grows to large positive values on one side and large negative values on the other side of the jet.

The observed values of $K$ are, inevitably, never as low as this extreme model would predict. This prompts a slightly more sophisticated model: instead of constant sea level to either side of the jet, add Gaussian noise to the sea level on either side. Thus, we consider the situation in which there is a sharp step in sea level across the jet, which meanders, but there is also random sea level variability producing a Gaussian PDF either side of the step (for simplicity, we choose the standard deviations of the Gaussians to be the same on each side of the step). The time series at any point is then a random variation about a certain mean at times when the point is to the north of the jet, and a random variation about a different mean at times when the point is to the south of the jet. The resulting PDF consists of the sum of two Gaussian PDFs, with centres separated by a distance $d$ (the size of the step). We measure $d$ in units of the Gaussian standard deviation, so that $d$ represents the ratio of step size to the size of the Gaussian noise. We now define $A$ as the integral of the Gaussian centered at the lower value (i.e. the probability of being on the low side of the step), and $B = 1 - A$ is the integral of the second Gaussian. Again, we will have $A = B = 0.5$ at the mean central position of the jet, with $A$ decreasing and $B$ increasing to one side of that point, and the converse on the other side.

At the jet centre, the sum of two Gaussians with small separation would produce a PDF with a flatter centre region than a single Gaussian alone (and hence kurtosis less than 3). For larger separation ($d > 2$) the PDF becomes bimodal, and kurtosis becomes even smaller, dropping towards 1 as the separation of the Gaussians becomes larger (when the noise to either side of the jet becomes a small fraction of the size of the step).

Fig. 7: Time series of relative vorticity (plotted as $zeta/f$) from two points with particularly low kurtosis, from (top) the Gulf Stream extension and (bottom) the Agulhas Return Current. Panels to the right show the corresponding PDF, calculated by kernel density estimation using a Gaussian kernel of standard deviation 0.2.

Fig. 7 shows the time series, and associated normalized PDF, for $zeta/f$ at two points of particularly low kurtosis in the Gulf Stream extension and the Agulhas Return Current. In both cases the vorticity clearly spends about half its time in the vicinity of one value, and half its time close to another value, with clusters about those two values overlapping to a degree, but clearly separate. This is reflected in the bimodal PDFs.

The model remains simplified, but this has the advantage that it can be solved analytically for skewness and kurtosis as a function of $A$ and $d$. The resulting relationships (derived in Appendix A) are

$$ K - 3 = frac{d^4P(1-6P)}{(1 + d^2P)^2} $$, (3) $$ S = frac{d^3P(A-B)}{(1+d^2P)^{3/2}} $$, (4) $$ S^2 = frac{d^6P^2(1-4P)}{(1+d^2P)^3} $$, (5)

remembering that $B = 1 - A$ and $P = AB$. At the mean position of the jet, $A = B = 1/2$ giving a skewness of zero, and reducing (3) to $K = 3 - 2/(1 + 4/d^2)^2$. A series of curves, each corresponding to changing $A$ at constant $d$ (i.e. representing how skewness and kurtosis would vary across a jet for a particular ratio of noise to step size) is shown as the thin black curves in Fig. 8. The different curves represent integer values of $d$ from 1 to 6 (innerÐouter), showing how the distance from the Gaussian values $S = 0$, $K = 3$ increases as the step size $d$ increases in comparison to the noise. To understand this plot, consider a jet with negative skewness to the south (for relative vorticity, an eastward jet in the southern hemisphere would show this pattern, for sea level it would be a northern hemisphere jet). As we move from south to north, we follow one of the curves in an anticlockwise direction. Far to the south, the noise dominates over the influence of the jet and the PDF is a single Gaussian ($P = 0$), producing $S = 0$, $K = 3$. As the influence of sea level from the other side of the jet becomes stronger ($P$ increases), skewness decreases and kurtosis increases, reaching a maximum at $P = 1/(d^2+12)$. Kurtosis then starts to decrease while skewness continues to decrease until it reaches a minimum at $P = 2/(d^2+12)$, where $K = 3 + 3S^2/2$. Kurtosis drops below 3 at $P = 1/6$, decreasing to its minimum at the centre of the jet where $P = 1/4$, $S = 0$. Then kurtosis starts to rise and the curve follows the positive skewness branch to the right and up, returning finally to $S = 0$, $K = 3$ far to the north.

Fig. 8: The relationship between skewness and kurtosis in the vicinity of a front. Thin curves represent the relationship resulting from the simple model described in the text. Each curve corresponds to a particular ratio d of step size to noise (values d .1Ð6 are plotted), with curves becoming larger as d increases (noise decreases). For relative vorticity, the curve is followed anticlockwise when passing from poleward to equatorward of an eastward jet. Crosses are values at points in the vicinity of the Agulhas Return Current. The thick line with diamonds is the trajectory described by statistics averaged at constant meridional distance from the centre of the current.

For small $d$ (small step size compared to noise), the trajectories remain close to $S = 0$, $K = 3$, extending further from this point as $d$ increases. In fact, allowing all values of $d$ produces a family of curves which Þll the entire allowed space $K geq S^2 + 1$. Only for $d > 2$ does the PDF become bimodal at the jet centre, at which point we would then have $K < 2.5$, although for more general distributions, the PDF can only be guaranteed to be bimodal if $K < S^2 + 1.512$ (Klaassen et al., 2000).

There are many ways in which this model is oversimplified. Noise is the same (and Gaussian) on either side of the jet, and the jet is still modelled as a step function. A slightly more realistic model has also been investigated numerically, modelling the jet as $a(t)tanh{(chi - chi_0(t))} + b(t)$, with Gaussian random time series for $a$, $b$, and $chi_0$. Results from this are not shown as the curves produced are very similar (though not identical) to those given by the simpler model, and do not seem to provide much greater insight. Two main criteria determine whether the trajectory forms a wide loop about the point $S = 0$, $K = 3$: the jet must meander by more than its width, and the sea level noise either side of the jet must be smaller than the sea level step across the jet.

We have illustrated the actual skewness-kurtosis relationship for the Agulhas return current in Fig. 8. The centre of the jet at each longitude between 23.3˚E and 50˚E was defined as the point of minimum kurtosis in a narrow latitude range south of Africa. Kurtosis and skewness were then regridded, shifting each longitude to the north or south so that it was centered on the jet centre. Skewness and kurtosis at all points from 10 grid points south to 10 grid points north of the centre (a total range of about 4.8 degrees) are plotted in the figure as grey crosses. Then skewness and kurtosis were averaged at constant northward distance from the jet centre, and the resulting trajectory plotted as the thick line in Fig. 8.

Although the relationship does not lie along a particular constant $d$ curve, it is clear that the trajectory forms a fairly broad loop around $S = 0$, $K = 3$, often outside the $d = 2$ contour, and with minimum $K$ below 2.5 consistent with the observed bimodal distribution, and with the interpretation as a meandering jet. The asymmetry is interesting, suggesting that noise is more significant in the Antarctic Circumpolar Current (ACC) to the south of the jet, and that the meanders are the more dominant signal on the northern side, in the Indian Ocean. Thus, although kurtosis falls below 3 only in a small region around the mean jet position, the statistics remain consistent with interpretation as the effect of a meandering jet to a distance of several degrees either side of the jet.

It will not have escaped the readerÕs attention that the discussion has been in terms of sea level, whereas most of the illustrations are of relative vorticity. While it is clear that a jet must be associated with a drop in sea level, the associated vorticity pattern is not so clear, except to note that, as velocity is a maximum at the jet centre (assuming an eastward jet), relative vorticity must be cyclonic on the poleward side and anticyclonic on the equatorward side of the jet, and small far from the jet. It is not inevitable that the strongest gradient in relative vorticity lies at the jet centre.

In one sense, it does not matter. The arguments above apply if there is a sharp, meandering front in any quantity, and the low kurtosis region can then be taken to represent the mean position of that front. However, there is a more concrete model which is rather appealing, which is clearest if we consider things in terms of potential vorticity (PV). One model of a jet, recently reviewed by Dritschel and McIntyre (2008, D&M hereafter, see also references therein), is as the result of a step in PV. In this paradigm, jets represent a barrier to the mixing of PV, so PV tends to become mixed to two different, constant values either side of a jet, with a sharp gradient across the jet. In turn, inverting the PV distribution shows that such a PV staircase actually implies the existence of a jet at the position of the PV step, so a positive feedback exists tending to maintain this configuration. As PV is materially conserved, time series of PV in the vicinity of the jet will be close to the ideal of switching between two constant values which produces the lowest possible kurtosis values for a given skewness.

Our time series is, of necessity, relative vorticity rather than PV. However, if we assume we are in a dynamical regime in which variations of PV are dominated by vorticity rather than thickness variations (or if we assume a constant relationship between thickness and relative vorticity changes), then the time series of relative vorticity at a fixed point (so that $f$ is constant) will be equivalent to time series of PV. In this case, we can see that the low kurtosis regions are consistent with the existence of a PV step at the jet centre. This prompts the question of whether kurtosis bears a relationship to the presence of mixing barriers, which will be considered in the next section.

First, though, we turn our attention to the many sharply-deÞned low kurtosis features which are clear in Figs. 5 and 6. The Gulf Stream, Kuroshio, and Agulhas Return Current are clear, but there are many other features, particularly in the Southern Ocean, which could also plausibly represent jets or frontal features. From our model, the centres of jets which meander strongly should be characterized by regions where the zero contour of skewness coincides with a value of kurtosis less than 3. In Fig. 9 we plot these contours in black, using $zeta/f$ skewness and kurtosis, and insisting on a degree of spatial continuity by using a 5 by 5 grid point smoothed version of the fields. The contours are superimposed on a map (top) of mean surface geostrophic ßow speed based on the Rio05 mean dynamic topography (Rio and Hernandez, 2004). It is immediately clear that the technique correctly identifies not only the major eastward-ßowing jets, but a number of weaker, mid-ocean jets in all ocean basins. There are, however, significant disagreements with the Rio05 data, most notably in the Southern Ocean. Here, agreement is generally good in the northern parts of the ACC, but is less good further south, especially in the eastern PaciÞc sector where a number of apparent jets are identified which are absent from Rio05.

However, small-scale features in the Rio05 data rely for their detection primarily on the surface drifter data, and these southern regions of the ACC are precisely where the density of drifter data drops off. It may be that there are really jets in these regions, but they are not resolved by the measurements available to produce the Rio05 map. In order to test this possibility, we have turned to sea surface temperature (SST) gradients. These are not always indicative of surface currents, but the fact that the ACC is to first-order an equivalent barotropic ßow (Killworth, 1992; Killworth and Hughes, 2002) means that temperature at any depth does tend to be a good proxy for the dynamic topography in this region (Sun and Watts, 2001). In order to calculate a mean SST gradient, we have used merged infrared and microwave satellite data from the Mersea Odyssea project, downloaded from ftp://ftp.ifremer.fr/ ifremer/medspiration/data/l4hrsstfnd/eurdac/glob/odyssea. The data are provided as daily fields at 0.1 degree resolution. Over the period 1st October 2007 to 12th May 2009, 451 fields were available. These were differentiated to produce temperature gradients and each field inspected by eye to check for obvious artifacts. On this basis, 30 daily fields were eliminated. The remaining 421 fields of SST gradient were averaged, and regridded onto the AVISO altimetry grid.

The resulting SST gradient plot is shown in the lower panel of Fig. 9, with our putative jet contours superimposed. It is clear that the Southern Ocean contours do indeed correspond with frontal features which were unresolved in the Rio05 data. In order to test how robust these comparisons are, we have repeated them with the Maximenko and Niiler (2005) dynamic topography, and with 7.25 years of SST data from the coarser-resolution AMSR-E microwave temperature measurements (Wentz and Meissner, 2004). The resulting pictures (not shown) are very similar, with a little smoothing of the SST gradient, and a little less smoothing of the geostrophic ßow speed (at the expense of more small-scale noise).

Fig. 9: Magnitude of the mean current (top) and mean sea surface temperature gradient (bottom), with black contours superimposed where the skewness/kurtosis relationship implies the mean centre of a jet or front to lie.

One final feature worth noting is that the identified path of the Kuroshio seems to be on the northern flank of the jet itself. The same is also true if sea level is used to identify this jet (not shown). This may be a result of the different averaging periods used for the mean current and the skewness and kurtosis; preliminary calculations using different subsets of the altimeter data produce more variation in results around the Kuroshio than elsewhere. Close inspection of the Gulf Stream also reveals some fine scale structure in the vorticity diagnostics which are not present in the sea level equivalents. It may be that, by combining the two analyses, more detailed information about the lateral structure of the jets can be extracted.

In summary, the pattern of a low kurtosis zone at the centre of a change in sign of skewness seems to be a reliable method of identifying meandering jets or frontal zones. When applied to time series of relative vorticity, this has allowed us to identify several jets previously unresolved by drifter data. In the next section, we will consider whether these low kurtosis zones also represent mixing barriers.

Does low kurtosis imply a mixing barrier?

As discussed above, the PV staircase model reviewed by D&M implies that jets occur at steps in PV, and that they act as mixing barriers. In fact D&M make the argument that mixing across a jet requires the presence of vortices with PV anomalies larger than the PV step which is associated with the jet. This sounds very similar to the criterion necessary for the formation of a low kurtosis region. If vortices are considered to be the ÔÔnoiseÓ on either side of the jet, then low kurtosis at the jet centre only occurs when the noise is smaller than the step, suggesting that low kurtosis should qualitatively be an indicator of a mixing barrier.

In order to test this suggestion, we turn to a set of idealized model experiments, initially run in order to investigate the effect of topography on jets and mixing in a zonal channel (Thompson, 2009). These are two-layer, quasigeostrophic simulations in a doubly-periodic channel, forced by an implicit zonal momentum flux which maintains a constant difference between the total zonal transports in the upper and lower layers. These model runs are not chosen because they are supposed to model particular parts of the real ocean, but because they explore a wide range of jet types and behaviors, including intermittent jets, meridionally-drifting jets, and topographically-steered jets, as well as steady, zonal jets. Another reason for the model choice is that the mixing diagnostics have already been performed, making it simple to compare mixing with kurtosis.

Four model experiments are considered. Experiment 1 is a ßat-bottomed channel, and produces four, steady, eastward jets, with weaker westward flow in-between. Experiment 32 has sinusoidal topography as a function of latitude only, and produces a sequence of jets which form in regions of strong background PV gradient, but drift north or south to regions of weak PV gradient, where they decay. Experiment 84 has topography which is sinusoidal in both latitude and longitude, producing a series of three strongly-steered eastward jets, with weaker westward ßow in between. Experiment 87 has similar topography, but only one strong eastward jet, which is more weakly steered and which intermittently breaks up in a bout of strong mixing. More details of these model runs are given in Thompson (2009). Mean and eddy kinetic energies vary by a factor of more than 40 between experiments, which clearly represent a wide range of conditions and kinds of jet.

Flow fields and PV kurtosis maps for the top layer are shown in Fig. 10 (in the case of experiment 32, this is a plot of zonally-aver aged ßow and kurtosis as a function of time and latitude, with kurtosis calculated using a moving window). It is clear that kurtosis is robustly low at the centre of the jets. For these experiments, the effective diffusivity $K_{eff}$ of upper layer PV was calculated as a function of time and equivalent latitude using the methodology of Shuckburgh and Haynes (2003). In order to compare across the range of experiments, it was found most effective to normalize the diffusivity by the square root of the domain-averaged EKE. The resulting normalized diffusivities would have dimensions of length, but the experiments were all non-dimensionalized using the Rossby radius $R_0$ as the length scale, and $R_0/U$ for time, where $U$ is the average imposed velocity difference between layer 1 and layer 2, so the normalized diffusivity should be multiplied by $R_0$ to give a dimensional value.

Fig. 10. Kurtosis of upper layer PV (shaded) and upper layer streamfunction (contours) for the four model experiments considered here. Plots are time averages, except for experiment 32, which is a zonal average as a function of time, the kurtosis being a running average over 20 time units.

The time-averaged values of $K_{eff}$, normalized by $sqrt{EKE}$, were compared with zonally-averaged values of PV kurtosis (in the case of experiments 84 and 87 the average was not zonal but along mean streamlines, because the mean flow is strongly-steered by topography). The resulting relationship between kurtosis and $K_{eff}/sqrt{EKE}$ is plotted for all experiments in Fig. 11. There is not just a qualitative relationship between the two, at low values of diffusivity there appears to be a strong, universal, linear relationship. As a guide (not any kind of formal Þt), the straight line $K_{eff}/sqrt{EKE} = (K - 1.3)/30$ has been drawn, which appears to fit the distribution fairly well.

Fig. 11: Scatter plot of surface layer normalized effective diffusivity against zonally-averaged kurtosis of PV (along stream-averaged for experiments 84 and 87) for the four numerical experiments. The solid line represents a suggested suitable linear relationship, but it is not a formal Þt. Values of the normalizing constant, non-dimensional domain-averaged $sqrt{EKE}$, are listed for each experiment.

For kurtosis values significantly higher than 3, the relationship becomes looser. There is no reason to expect a relationship at these high $K$ values, so this is unsurprising. The relationship is also significantly worse for experiment 32, which again is not surprising as the value at each latitude is an average over times when there is low mixing and times when there is high mixing at that latitude. The surprise is that it works as well as it does, and falls on roughly the same range of values as the other experiments. Fig. 12: The same values of normalized effective diffusivity from Fig. 11, plotted as a function of effective latitudinal coordinate measured in Rossby radii. The dashed line is the value which would be predicted from kurtosis of PV, using the linear relationship suggested in Fig. 11.

Fig. 12 uses the linear scaling proposed in Fig. 11 to plot normalized $K_{eff}$ as a function of latitude for each experiment, together with the value which would be predicted from the kurtosis (dashed lines). This shows how well kurtosis captures the structure as well as values at the minima of diffusivity. There is a tendency to underestimate diffusivity where jets are drifting (experiment 32) or intermittent (experiment 87), in which case our simple model is no longer adequate, but even here there is a reasonable correspondence. Much better fits can be found for individual experiments by adding quadratic terms, but that is at the expense of universality, since the quadratic terms differ strongly between experiments.

Fig. 12: The same values of normalized effective diffusivity from Fig. 11, plotted as a function of effective latitudinal coordinate measured in Rossby radii. The dashed line is the value which would be predicted from kurtosis of PV, using the linear relationship suggested in Fig. 11.

Conclusions and discussion

Meandering jets in the ocean are associated with a characteristic pattern of skewness and kurtosis of sea level and, more clearly, relative vorticity. The mean position of the jet or front lies along the zero contour of skewness, which is also a region of low (less than 3) kurtosis. This can be explained using a model in which the jet or front represents a sharp step in sea level or relative vorticity, and meanders over a distance wider than the width of the step. We find little evidence that more useful data can be extracted from the spatial variation of PDFs beyond what is apparent from skewness and kurtosis. Using this relationship we have identified several Southern Ocean jets which were previously unresolved by drifter-based climatologies.

The fact that these features are seen more clearly in relative vorticity than in sea level suggests the use of a model of a jet as a step in PV, between two regions of well-mixed PV, as discussed recently by D&M. In this model, when the jump in PV across the step is larger than the eddy variations to either side, the step acts as a mixing barrier. These are precisely the conditions in which time series of PV would exhibit low kurtosis at the mean position of the jet centre. That led us to ask whether low kurtosis regions were indicative of the presence of a mixing barrier.

Across a wide variety of idealized model experiments, we find a strong, universal, linear relationship between PV kurtosis and effective diffusivity normalized by the square root of domain-averaged EKE. The relationship is particularly strong at low values of kurtosis, which represent the strongest mixing barriers. This suggests that the regions of low kurtosis from the relative vorticity time series may indeed represent mixing barriers in the ocean.

There are a couple of difficulties with linking the model experiments and our interpretation to the data. The data only tell us about relative vorticity, not PV. Variations in $f$ are not a problem, as the time series are at constant latitude, but if thickness variations are not either small or simply related to relative vorticity variations, then the time series we have may not be representative of PV variations. The second problem is that the relationship we find is between PV kurtosis and a diffusivity normalized by the square root of domain-averaged EKE. In the real ocean, what would be the appropriate region over which to average? Using local EKE in the model diagnostics results in significantly more scatter in the relationship, so some degree of spatial averaging is clearly needed, perhaps over an eddy Rhines scale? These difficulties mean that it would be premature to produce quantitative estimates of real ocean diffusivities using this method: a study is needed using more realistic ocean model geometry to determine how best to apply the relationship we find.

Kurtosis does, however, have two advantages over effective diffusivity: it is very simple to calculate, and it is a local quantity, meaning that maps are easily produced, unlike effective diffusivity which only varies in one spatial dimension (although it is less obvious how to produce time series of kurtosis than of diffusivity). This means that it may be useful for quick identifications of mixing barriers in a number of contexts. For example, it would be simple to extend the present analysis to subsurface regions in an ocean model, permitting the investigation of the three-dimensional geometry of mixing barriers. Similarly, the ideas should apply equally well to atmospheric data. As White (1980) showed, a similar pattern (a zero skewness contour between two regions of significant skewness of opposite signs, coinciding with a region of low kurtosis) occurs in maps derived from 500 mbar geopotential height, though not at 1000 mbar. As with sea level, it may be that this picture would be sharpened by considering relative vorticity or PV.

There is a good case to make skewness and kurtosis (whether it is of sea level, relative vorticity, or PV) standard diagnostics to be used to assess the realism of eddying ocean models, as well as to understand the dynamics of jets and mixing.

Appendix A. Derivation of statistics for a PDF which is the sum of two Gaussians

If we write the Gaussian distribution with standard deviation $sigma$ as

$$ G(chi, sigma) = frac{1}{sigmasqrt{2pi}}expleft({-chi^2/2sigma^2}right) $$, (6)

which has been normalized so that its integral is 1, then integration by parts shows that

$$ int_{-infty}^{infty}{chi^2 , G(chi, sigma) dchi} = sigma^2 $$, (7)

$$ int_{-infty}^{infty}{chi^4 , G(chi, sigma) dchi} = 3sigma^4 $$, (8)

and odd moments of $G$ are zero by symmetry.

The PDF which consists of two Gaussians separated by $d$ standard deviations is

$$ p = A G(chi + Bsigma d, sigma) + B G(chi - A sigma d, sigma) $$, (9)

where $B = 1 - A$, and the origin has been chosen so that the mean value of the PDF is zero (to construct this, note that the distances of the Gaussian centres from the origin must be inversely proportional to their integrals $A$ and $B$, and that the distance between the Gaussians is deÞned to be $sigma d$). The second moment of this distribution, which we will set to 1, is then given by

$$ int_{-infty}^{infty}{chi^2 , p dchi} = A int_{-infty}^{infty}{(chi - B sigma d)^2 G(chi, sigma)} dchi + B int_{-infty}^{infty}{(chi + A sigma d)^2 G(chi, sigma)} dchi $$. (10)

Expanding the squared terms, odd terms in $chi$ integrate to zero, and even terms were evaluated above, so this gives, after some gathering of terms,

$$ 1 = sigma^2(1 + d^2P) $$, (11)

where $P = AB$, or

$sigma = (1 + d^2P)^{-1/2}$. (12)

The same procedure can be applied to the third and fourth moments, leading to

$$ S = sigma^3 d^3 P(A - B) $$, (13)

and

$$ K = 3sigma^4 + d^2 sigma^4 P[6 + d^2(1 - 3P)] $$. (14)

Substituting for $sigma$ and simplifying then leads to (3) and (4).