The statistical nature of turbulent barotropic ocean jets

Jets are an important element of the global ocean circulation. Since these jets are turbulent, it is important that they are characterized using a statistical framework. A high resolution barotropic channel ocean model is used to study jet statistics over a wide range of forcing and dissipation parameters. The ﬁrst four moments of the potential vorticity distribution on contours of time-averaged streamfunction are considered: mean, standard deviation, skewness and kurtosis. A self-similar response to forcing is found in the mean and standard deviation for eastward barotropic jets which exhibit strong mixing barriers; this self-similarity is related to the global potential enstrophy of the ﬂow. The skewness and kurtosis give a behaviour which is characteristic of mixing barriers, revealing a bi/trimodal statistical distribution of potential vorticity with homogenized potential vorticity on each side of the barrier. The mixing barrier can be described by a simple statistical model. This behaviour is shown to be lost in westward jets due to an asymmetry in the formation of zonal mixing barriers. Moreover, when the statistical analysis is performed on eastward jets in a streamfunction following frame of reference, the distribution becomes monomodal. In this way we can distinguish between the statistics due to wave-like meandering of the jet and the statistics due to the more diffusive eddies. The statistical signature of mixing barriers can be seen in more realistic representations of the Southern Ocean and is shown to be an useful diagnostic tool for identifying strong jets on isopycnal surfaces. The statistical consequences of the presence, and absence, of mixing barriers are likely to be valuable for the development of stochastic representations of eddies and their dynamics in ocean models.


Introduction
Ocean jets such as the Gulf Stream, Kuroshio Extension, as well as the jets embedded within the Antarctic Circumpolar Current, represent some of the most turbulent regions of the ocean (e.g., Stammer and Wunsch, 1999 ). Such jets have an important influence on the large scale circulation of the ocean and atmosphere. For example, Scaife et al. (2011) highlight that model representation of winter atmospheric blocking is greatly improved when the horizontal resolution of the North Atlantic ocean is increased, thereby improving the representation of the Gulf Stream. From a hydrographic survey, Bower et al. (1985) determined that water masses are mixed across the Gulf Stream at depth, but that the Gulf Stream acts as a barrier towards the surface where the jet is strongest, inhibiting exchange of tracers between the subtropical and subpolar gyre. More recent studies ( Ferrari and Nikurashin, 2010;Klocker et al., 2012 ) show the suppression of mixing, and consequently the transport of heat and salt, across jets within the Antarctic Circumpolar Current. This suppression of mixing can be attributed to the tendency of strong jets to form mixing barriers, regions of high potential vorticity gradient which eddies have difficulty penetrating.
The tendency for strong jets to form mixing barriers is, in fact, a consequence of the turbulence which sharpens the jets and maintain strong potential vorticity gradients by providing up-gradient fluxes of momentum ( Starr, 1968 ). Evidence for this jet sharpening behaviour has a long history in baroclinic (e.g. Fultz et al., 1959 ) and barotropic (e.g. Sommeria et al., 1989 ) jets which is comprehensively reviewed in Dritschel and McIntyre (2008) . The highly anisotropic picture we have formed from these studies is that of a potential vorticity step corresponding to the jet itself, sharpened by the eddies, with regions of near homogenized potential vorticity to either side of the jet within which the eddies act diffusively.
The dynamical description of ocean jets is complemented by attempts to derive statistical models of turbulent geophysical flow ( Esler, 2008 ). This is motivated in part due to the stochastic nature of turbulent flows which necessarily require more than a purely deterministic description. Moreover, with increasing interest in the application of stochastic parameterization to ocean modelling (e.g., Right : examples of leptokurtic ( red ), kurtosis greater than three, and platykurtic ( blue ), kurtosis less than three, distributions compared with a normal distribution with a kurtosis of three ( black ). All distributions shown have a mean of 0 and a standard deviation of 1. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Berloff, 2005;Porta Mana and Zanna, 2014;Grooms et al., 2015;Zanna et al., 2017 ) it is important to understand the underlying statistics of geophysical flows and the fundamental organizing principles behind them. Many statistical mechanics theories of turbulent flow are derived in the absence in of forcing and dissipation: the so-called Robert-Sommeria-Miller equilibrium statistical mechanics ( Bretherton and Haidvogel, 1976;Salmon et al., 1976;Robert, 1990;Miller, 1990;Robert, 1991;Robert and Sommeria, 1991 ). Equilibrium statistical mechanics theories of geostrophic turbulence predict a functional dependence of the mean potential vorticity on the mean streamfunction (e.g., Bouchet and Venaille, 2012 ). However, the ocean lies in a forced-dissipative regime for which a non-equilibrium statistical mechanics theory is required.
The goal of this study is to understand the statistics of forceddissipative ocean jets and how these are influenced by the dynamics of a mixing barrier. We explore the sensitivity of the jet statistics to forcing and dissipation and show that self-similarity exists across a wide range of parameters, describing the characteristic statistical signature of a mixing barrier.
To address this goal, we examine the mean state of the flow and its variability in the presence of forcing and dissipation using mean, standard deviation and higher order moments such as skewness and kurtosis. We will also compare and contrast these statistics to the predictions of equilibrium statistical mechanics. The energy and enstrophy conserving theory considered in this study is that presented in Jung et al. (2006) . Jung et al. (2006) differs from other tests of equilibrium statistical mechanics (e.g. Qi and Marston, 2014;Dritschel et al., 2015 ) in that it derives a spatially coarse-grained version of the theory in order to test the theory in a barotropic rotating annulus laboratory experiment. When only energy and potential enstrophy are conserved, the relationship between mean potential vorticity and the mean streamfunction is linear while the standard deviation is constant with respect to streamfunction according to equilibrium statistical mechanics. These predictions are derived in Jung et al. (2006) and are consistent with other interpretations of the theory (e.g. Salmon et al., 1976 ). Understanding these relationships in a forced-dissipative system, however, remains an open problem which we will discuss during the course of this study.
In addition to understanding how the mean and standard deviation relations change in the presence of forcing and dissipation, it is important to consider higher order moments of the flow, such as skewness and kurtosis, which exist only when the potential vorticity distribution is non-Gaussian. As illustrated schematically in Fig. 1 , different values of skewness and kurtosis of the probability distribution characterize asymmetry and intermittency (extreme events) in turbulent flow. These moments have been shown to be important in describing the statistics of passive tracers (e.g., Bourlioux and Majda, 2002 ). The character and nature of skewness and kurtosis in turbulent geophysical flows have been considered by Thompson and Demirov (2006) and Hughes et al. (2010) . In particular, Hughes et al. (2010) establish a mechanistic link between skewness and kurtosis and the large scale dynamics of jets, allowing these authors to use these moments to identify mixing barriers. In a highly non-linear dynamical system such as the ocean, we might expect extreme events in active tracers such as potential vorticity to have an impact on the mean large-scale circulation and to be sensitive to forcing and dissipation. In this manuscript, we study the potential vorticity statistics of a barotropic jet in a reentrant channel, varying the strength of forcing, the linear drag coefficient, and the direction of forcing. These experiments are used to pursue the following aims: • To understand the dependence of the statistics of a turbulent flow on forcing and dissipation in the presence and absence of a mixing barrier. • To establish a statistical picture of mixing barriers with a simple statistical model which will inform the stochastic parameterization of eddy-mean flow interaction. • To distinguish between the statistics of eddies and the statistics of waves.
The results are restricted here to barotropic flows, an understanding of which is a prerequisite for understanding more complex flows with baroclinicity.
The structure of the paper is as follows. In Section 2 , the model and experiments are presented. In Section 3 , the dynamics of the west/eastward jet asymmetry is discussed. In Section 4 , the mean and standard deviation of potential vorticity are calculated as a function of time-averaged streamfunction in order to contrast with the equilibria of the ideal dynamics, and a self-similar response to forcing and dissipation is discussed. In Section 5 , the higher order moments of skewness and kurtosis are similarly considered. A clear link is made between the presence of a meandering mixing barrier and a characteristic skewness/kurtosis behaviour. In Section 6 , the importance of the meandering mixing barrier in determining the statistics of the flow is further considered, and the statistics calculated on instantaneous contours of streamfunction in order to follow the large scale meanderings of the mixing barrier in a quasi-Lagrangian sense. In Section 7 , a further discussion of: joint potential vorticity and streamfunction statistics; the importance of mixed distributions; implications for stochastic parameterization; and application to Southern Ocean jets using an observation constrained primitive equation model. Finally, Section 8 contains some concluding remarks.

Model and dynamics
In this study we focus on a idealized barotropic model of forced-dissipative turbulent flow in a channel. Although ocean jets are baroclinic in nature, barotropic instability plays an important role in the large scale dynamics of the meandering jet ( Waterman et al., 2011 ). In this way, by learning about the statistical nature of turbulent jets in the idealized system considered here, we can build expectations of the behaviour of more complex systems. Recent studies on the equilibration of freely decaying turbulence are Esler (2008) and Harnik et al. (2014) which consider the baroclinic and barotropic channel problem respectively. In both studies the formation and statistics of mixing barriers are important in determining the equilibrated state. Here we consider the forceddissipative statistics. In addition, the use of a barotropic model allows us to draw a direct comparison of Jung et al. (2006) who found good agreement between a barotropic rotating annulus experiment and an (energy-enstrophy) equilibrium mechanics theory.
The non-dimensional, prognostic equation to be solved is where ψ is the streamfunction, is the horizontal Jacobian operator, τ ( y ) is the zonal wind stress which varies only with meridional coordinate y, r is the linear drag coefficient and ν h is the hyper-viscosity. Linear drag is included to remove momentum at large scales and the hyper-viscosity is included in order to preferentially damp potential enstrophy at small scales. The potential vorticity is given by where ∇ 2 ψ is the relative vorticity and βy is the planetary vorticity. Here we adopt a β-plane approximation in which the planetary vorticity varies linearly with y . In this study we refer to q as the potential vorticity, and to the integral of q 2 as the potential enstrophy, in order to facilitate comparison with studies that incorporate baroclinicity; however, in the barotropic model the potential vorticity is identical to the absolute vorticity. The numerical code used is a modified version of that developed for the freely decaying turbulence problem in Esler (2008) and, originally, Esler and Haynes (1999) . The model is solved using a quasi-spectral method, spectral in the zonal direction making use of the periodicity. In the meridional direction a finite difference method is used. The first obvious difference is that the original code solved for two-layer quasi-geostrophic dynamics whereas the modified code solves for a barotropic flow with rigid lid. Secondly, a wind stress forcing has been added which required an alternative boundary condition to be applied to ensure the solution is consistent with momentum conservation. The domain geometry is given in Fig. 2 . The boundary conditions follow those specified in McWilliams (1977) requiring an integral momentum balance to be applied at each time-step. The no normal flow condition is written where the subscripts refer to the northern and southern boundaries respectively. To calculate the hyper-viscosity term we need two boundary conditions. To ensure no boundary stress we use the free slip conditions, where n = 1 , 2 . Finally, the momentum conserving condition is written where ( t ) is spatially constant along the north and south boundaries and evolves temporally through the integral momentum balance, Two parameters are varied across a wide range of experiments, labelled A through to U. In the first group of experiments, A-G, the strength of forcing, τ 0 , is changed holding the linear drag coefficient, r , constant. The meridional variation of the wind stress is kept the same for all experiments, where δ is a constant controlling the width of the wind forcing. In the second group of experiments, H-N, the linear drag coefficient, r , is changed while holding the forcing strength, τ 0 , constant. Finally the third group of experiments is a repeat of the first but with the direction of the forcing reversed. The parameters for all three experiments are given in Table 1 with the parameters which are held constant for all experiments are given in Table 2 .

Linear stability and west/eastward asymmetry
During this study we will refer to jets as any peaked zonal mean flow whereas the term mixing barrier will refer to any region of high potential vorticity gradient which consequently prohibits transport across it. A major dynamical difference is seen between the jets in simulations A-G (with varying forcing strength) and the jets in simulations O-U, with reversed forcing. This difference is evident in the snapshots of potential vorticity for simulations D and R, shown in Fig. 3 for illustrative purposes as the flow is typical for the respective experiments: A-G and O-U. The main difference to observe is the clear loss of the narrow region of high cross flow potential vorticity gradient along with the concomitant large scale Rossby wave. This is a persistent structure in the eastward jet case for all forcing strengths but is not at all apparent in the westward propagating jet. It would not be correct to say that mixing barriers do not form at all in the reversed forcing experiments, this can be seen in a snapshot from simulation U in Fig. 4 . While there is no barrier prohibiting mixing across the channel, we do see mixing barriers which wrap themselves around two large coherent eddies, partially isolating potential vorticity within the two eddies. The topologically distinct nature of these mixing barriers changes the barriers' influence on the statistics as a function on time-averaged streamfunction, which is zonally symmetric. A schematic representing the various flow regimes simulated in this study is given in Fig. 5 . This drastic difference between the two experiments implies that the dynamics are not symmetric with respect to a change in the sign of forcing. It can be seen from Eqs. (1) -(3) that changing the sign of forcing is equivalent to changing the sign of β combined with a reflection in the meridional ( y ) axis. Making such a transformation leads to a drastic difference in the resulting statis-tically steady state flow. The presence and absence of this mixing barrier will form a large part of the discussion in the following sections.
The Rayleigh-Kuo inflection-point criterion ( Rayleigh, 1880;Kuo, 1949 ) states that a necessary condition for barotropic instability of a zonal flow is that the meridional potential vorticity gradient changes sign, i.e., somewhere in the domain where U ( y ) is the zonally symmetric velocity profile of the mean jet. While (9) is a necessary global condition for instability, in practice, instability often occurs in the locality that (9) is met (e.g., Waterman and Jayne, 2012;Waterman and Lilly, 2015 ). During spin up, the zonally symmetric forcing profile strengthens a jet with the same meridional profile until (9) is satisfied; in the case of the westward jet, (9) is first satisfied at the centre of the channel ( y = 0 ); whereas in the case of the eastward jet, (9) is first satisfied on jet flanks. It is in these very same regions that the nonlinear eddies emerge and erode the background potential vorticity gradient in Fig. 4 . Thus, while the linear instability condition is not valid in the fully nonlinear, turbulent regime, it is nonetheless successful in predicting the initial regions of eddy stirring and subsequent emergence of mixing barriers in the adjacent regions where the background potential vorticity gradient and Rossby restoring force is enhanced.
In following sections, the statistics of the eastward and westward jets are compared and contrasted in order to isolate effects associated with strong mixing barriers.

Scaling of the mean and standard deviation of potential vorticity
In this section we consider the relationship between the mean potential vorticity, q , and the time-averaged streamfunction, ψ , Fig. 3. Snapshots of potential vorticity from the statistically steady state for simulations D and R. Simulation D, eastward wind stress, has strong cross-channel potential vorticity gradients forming a mixing barrier which meanders due to persistent travelling Rossby waves. Simulation R, westward wind stress, does not show the same behaviour rather exhibits mixing across the channel.  as well as the relationship between the standard deviation of potential vorticity and the time-averaged streamfunction. The angled brackets, , imply an average over all points in time and space coincident with a contour of time-averaged streamfunction. For the numerical simulations presented here it is important to note that the time-averaged streamfunction is zonally symmetric since the simulations are simulation for a sufficiently long time in steady state. Therefore there is a one to one mapping between meridional distance and time-averaged streamfunction. In We see a characteristic functional dependence of both the mean and the standard deviation on the time-averaged streamfunction.
There is a stark contrast between these relations and those predicted by the equilibrium energy-potential enstrophy statistical mechanics: mean linear, standard deviation constant ( Salmon et al., 1976;Jung et al., 2006 ). However, the mean potential vorticity dependence on streamfunction bears resemblance to the tanhlike relations ubiquitous in the generalized equilibrium Robert-Sommeria-Miller theory ( Venaille and Bouchet, 2011 ). This implies that these statistics require more than knowledge of the energy and potential enstrophy to be explained. A detailed consideration for the applicability of statistical mechanics theory to forceddissipative flow is left for future studies. Here we highlight the generic nature of these statistics and their behaviour over a wide range of forcing and dissipation. Fig. 6 c and d shows that the mean and standard deviation of potential vorticity display a self-similar dependence on the normalized streamfunction when both the q − ψ relation and the SD (q ) − ψ relation are rescaled by their 2-norm. The red dashed lines in Fig. 6 c and d gives the error in this self-similarity which is found to be small. Interestingly it is possible to relate the nor- Fig. 6. a) mean potential vorticity against normalized time-averaged streamfunction. b) standard deviation of potential vorticity against normalized time-averaged streamfunction. The streamfunction is normalized by it's steady state boundary value, , which is the integrated zonal momentum. Each plot shows the results from seven experiments each corresponding to a doubling of the forcing strength from the previous. The data from panels a and b are rescaled with respect to their norm to give; c) the average across all seven experiments of the rescaled mean potential vorticity as the solid-black curve; d) the average across all seven experiments of the rescaled standard deviation of the potential vorticity as the solid-black curve. In both c and d the dashed-red curves give the standard error. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) malization factors for q and SD( q ) to the potential enstrophy of the flow. Assuming zonal symmetry we can define a mean flow potential enstrophy as and an eddy potential enstrophy as and L x is the zonal extent of the domain. For the mean, the normalization factor is given as while for the standard deviation the normalization factor is given as We find that the q − ψ relation and the SD (q ) − ψ relations are self-similar over a wide range of eastward forcing strengths ( Fig. 6 ), but not when the magnitude of the dissipation is varied, or the sign of the forcing is changed, causing the jet to flow in the opposite direction ( Fig. 7 ). The perturbed dissipation experiments, if anything, show limited response in the mean while there is a scaled, although less clean, response in the standard deviation of the potential vorticity. We suggest that this difference is related to the fact that scale dependence of the forcing is unchanged when its strength is varied, it always acts at the first zonal wavenumber. On the other hand, the linear drag term interacts strongly with the flow over a different range of scales as the linear drag coefficient is varied. Investigating this effect would require additional experiments in which the width of the wind stress profile is varied as well as the structure of the wind forcing; this is not explored further here. An interesting feature of the experiments in which the dissipation is varied ( Fig. 7 a and b) is the asymmetry in the SD (q ) − ψ relation about ψ = 0 for the low dissipation simulations, in particular H and I. We can see a clue to the resolution of this puzzle in Fig. 8 which shows a snapshot of steady state potential vorticity for simulation H, the simulation with the lowest linear drag coefficient. We find that a large coherent eddy has emerged in bottom half of the domain. A similar eddy can also be seen in simulations I and J, see Fig. 5 ; in the case of J the eddy appears in the upper part of the domain. It is unclear why this eddy forms when the friction is decreased. When the direction of the forcing is reversed, the mean and standard deviation dependence on the time-averaged streamfunction is drastically altered and again there is some self-similarity, although this is much less clear that for the eastward jets. The main physical difference in the reversed forcing experiments is the loss of the clear mixing barrier formed through Rossby wave dynamics. When the forcing is reversed the instabilities are strong enough to ensure that a persistent Rossby wave is broken up and this changes the nature of the flow in a significant manner.  In summary, it is found for eastward forcing that the mean and standard deviation of the potential vorticity scales with forcing and that the dependence of the mean and standard deviation on the time-averaged streamfunction is self-similar over a large range of forcing strengths. The q − ψ relation is self-similar when the streamfunction is scaled by the total through channel momentum and the potential vorticity is scaled by the square root of the mean field potential enstrophy. The SD (q ) − ψ relation is selfsimilar when the streamfunction is scaled by the total through channel momentum and the potential vorticity is scaled by the square root of the eddy potential enstrophy. Thus, for this experimental set-up it is possible to empirically deduce the mean and standard deviation given that we know the potential enstrophy behaviour of the flow. It is inconclusive whether the mean scales with dissipation but there is clear response in the standard deviation to changed dissipation, although self-similarity is not seen. We see a changed behaviour in the reversed forcing experiment and leave a full description of this flow regime for future work. In the following sections we concentrate on describing the statistical phenomenology for the forcing experiments, A-G, referring to the other experiments when pertinent to the argument.

Skewness, kurtosis and jet meandering
Gaussian statistics can be characterized using only knowledge of the mean and standard deviation. When the statistics are non-Gaussian it becomes important to consider higher order statistics. We will here consider the 3rd and 4th order centred and standardized moments, the skewness and kurtosis, 1 given respectively as Fig. 9. a) skewness of potential vorticity as a function of time-averaged streamfunction. The line S = 0 is shown and represents statistics which is symmetric about its mean. b) kurtosis of potential vorticity as a function of time-averaged streamfunction. The lines K = 3 and K = 1 are shown. A kurtosis of 3 corresponds to a Gaussian distribution, a kurtosis greater than this is referred to as leptokurtic and while a kurtosis of greater than three is referred to as platykurtic. A kurtosis of 1 is the smallest possible kurtosis and corresponds to an unbiased coin toss. In c and d the skewness and the kurtosis are given respectively from the bimodal distribution, (17) , and the trimodal distribution, The skewness is a measure of how asymmetric or lopsided a statistical distribution is about its mean while the kurtosis is a measure of 'peakiness'. A high kurtosis implies a distribution which has a sharp peak with fat tails, a low kurtosis implies small tails but a flat peak. 2 Fig. 1 illustrates these moments schematically. Plots of the skewness and kurtosis against normalized time-averaged streamfunction for the simulations A-G are given in Fig. 9 a and b. Contrasted with the scaling behaviour of the mean and standard deviation, for central values of the normalized time-averaged streamfunction, the skewness and kurtosis display invariance with respect to changed forcing. Not only is the skewness and kurtosis behaviour in the jet core behaviour fixed with respect to forcing strength, but this behaviour can interpreted employing the results of Hughes et al. (2010) and showing how they are influenced by forcing strength, dissipation of reversed jet direction. Hughes et al. (2010) show that a change in sign of skewness coincident with low kurtosis values is indicative of a mixing barrier. This is shown both in relative vorticity anomaly data inferred from altimeter observations as well as in a two-layer quasi-geostrophic model in a doubly periodic channel. In this study, the same skewness/kurtosis behaviour is seen in a barotropic model in a singly-periodic channel.
Central to the argument of Hughes et al. (2010) is the idea that a jet which acts as a mixing barrier leads to a step in potential 2 This picture relies on having an uni-modal distribution with one peak. vorticity with homogenized values of potential vorticity either side of the barrier. This potential vorticity staircase paradigm for jets is reviewed in Dritschel and McIntyre (2008) and paints a picture (along with the references therein) in which mixing or homogenization of potential vorticity either side of a jump in the value of potential vorticity acts to reinforce the mixing barrier and to amplify the jump in potential vorticity. This feedback leads us a picture where an observer standing at the centre of the jet, defined by either ψ = 0 or y = 0 , would see the value of potential vorticity flip from negative to positive as the mixing barrier meanders past due to non-linearities or Rossby waves.
With the potential vorticity step picture in mind, we can explain the skewness/kurtosis behaviour shown in Fig. 9 . A kurtosis of K = 1 is described as sub-Gaussian. Kurtoses of greater than three are described as leptokurtic and kurtoses of less than three are described as platykurtic. In Fig. 9 , the skewness plot is marked with the line S = 0 . The kurtosis plot is marked with the K = 3 line corresponding to a Gaussian distribution and the K = 1 line which corresponds to an unbiased coin toss or the unbiased Bernoulli process, a stochastic process where the random variable can take one of two values. We see a distinct change in the value of skewness for a time-averaged streamfunction of around zero. The skewness behaviour is consistent across all forcing strengths for a range of −0 . 4 < ψ / < 0 . 4 , which we refer to as the core of the jet.
For larger values of normalized time-averaged streamfunction the skewness becomes more varied across the different simulations. Coincident with the described behaviour of skewness we see similar consistency in the kurtosis for the range −0 . 4 < ψ / < 0 . 4 . The kurtosis in the core of the jet is predominately platykurtic and we also see the kurtosis becoming close to sub-Gaussian. This skewness/kurtosis behaviour is consistent with the findings of Hughes et al. (2010) and can be interpreted as due to a flip in the value of potential vorticity as the mixing barrier meanders, through the action of passing Rossby waves.
In fact, we can go further than the above discussion and consider a simple analytical model for the skewness and kurtosis used by Hughes et al. (2010) . This model describes the statistics as Gaussian either side of a step in the value of potential vorticity, with the meandering of the jet represented as a Bernoulli process. Thus, the centred probability distribution for the potential vorticity is given by the sum of two normal distributions scaled by factors, a and 1 − a and separated by d times their standard deviation, σ . a is a decreasing function of meridional distance to represent the decreasing probability of observing positive potential vorticity as the jet is crossed. The probability distribution is is the normal distribution as a function of potential vorticity, q , with a mean, μ, and standard deviation, σ . By making the mean of each of the Gaussian peaks depend on a we ensure the total distribution has zero mean, giving us the centred statistics we need to calculate the skewness and kurtosis. The skewness and kurtosis corresponding to this distribution is shown in Fig. 9 c and d respectively. These analytically derived curves bear a striking resemblance to the observed behaviour as a function of streamfunction shown in Fig. 9 a and b (and also for the corresponding curves in Fig. 10 for the dissipation experiment only). This simple model is able to qualitatively reproduce the skewness/kurtosis behaviour extending away from the centre of the jet explaining the high peaks in kurtosis at large values of streamfunction -as one of peaks gets very small the kurtosis it appears as a very large skewed tail of the larger peak leading to a high value of kurtosis. However, as the height of the now small peak goes to zero, the kurtosis must go to three -the kurtosis of a normal distribution. Due to the close similarity between the analytical curves and the model diagnosed curves, it tempting to forget that the former is plotted as a function of normalized time-averaged streamfunction, ψ / , while the later is a function of the weighting factor, a . Thus the resemblance suggests a close relationship between a and ψ / , both monotonic functions of meridional distance.
One matter of concern lies in the fact that when we plot the distribution of potential vorticity for a given time-averaged streamfunction the statistics is trimodal, Fig. 11 , unlike the bimodal distribution, (17) . To take account of this, the calculation of Hughes et al. (2010) is modified by a central Gaussian peak of height weight, c , giving a probability distribution of The results of this calculation is also shown in Fig. 9 . The trimodal nature of the statistics observed in some of the experiments will be discussed further in Sections 6 and 7 . However it is clear from The scale for the distribution above is omitted as it can be normalized away and we are only interested in the shape of the distribution at present. Fig. 9 that the addition of a central peak does not change the qualitative story.
In a barotropic channel, we have observed a change in sign of skewness across a jet, this, combined with a low kurtosis, is indicative of a mixing barrier. A highly simplified model describing the statistics of the jet meandering as a Bernoulli distribution with superimposed small-scale Gaussian noise gives surprisingly good qualitative agreement with the observed skewness and kurtosis across the jet over a wide range of forcing and dissipation. Fig. 10 shows how this behaviour changes in the dissipation and reversed forcing experiments. For the dissipation experiments, Fig. 10 a and b, we see a similar picture to the forcing experiments for all but the two lowest dissipation parameters. These outliers are linked to the clear asymmetry in the standard deviation shown in Fig. 7 b. For the other dissipation parameters we see a crossing of the S = 0 line coincident with platykurtic kurtoses close to K = 1 . The real difference is seen in the reversed jet experiment ( Fig. 10 c and d). Here the skewness behaviour is in the opposite sense, now going from positive to negative, and the values differ greatly between the different (reversed) forcing strengths. In the core of the jet the kurtosis is leptokurtic for all simulations and does not display the near sub-Gaussian behaviour in the other two experiments.
The drastic change of the skewness/kurtosis behaviour can be understood when we realize that in reversing the direction of forcing, resulting in a east to west rather than a west to east jet, we have changed the jet dynamics, as discussed in Section 2 , causing the mixing barrier to be broken up. This allows the potential vorticity to be mixed in the north-south direction, resulting in a jet without a mixing barrier, destroying the coin-toss like picture and the corresponding signal in the skewness/kurtosis behaviour. When the mixing barrier is present the skewness/kurtosis behaviour can be explained by a mixed distribution; this will be discussed further in Section 7 . We have shown here that in a barotropic re-entrant channel that the presence of a mixing barrier results in a characteristic skewness/kurtosis across jet behaviour. When the mixing barrier is removed, the skewness/kurtosis behaviour does not persist.

Streamfunction following statistics
In the previous section we saw that the skewness/kurtosis behaviour can be partially explained through a simple statistical description in the presence of a mixing barrier and that the observed skewness/kurtosis behaviour is indeed characteristic of a mixing barrier. In this section, we consider the dependence of the statistics on the coordinate system in which the statistics calculated. We consider three frames of reference, leading to three methods of time-averaging: • Eulerian time-averaging, in which the statistics are calculated along the zonally-symmetric, time-mean streamfunction contours. • Galilean transformed time-averaging, in which the statistics are calculated along the non-zonally symmetric time-mean streamfunction contours in a frame of reference moving at the phase speed of the dominant Rossby wave. • Streamfunction following time-averaging, in which the statistics are calculated along the instantaneous streamfunction contours.
In the interests of brevity, we present results only for the experiments with the changed forcing strength, A-G.
We start by considering in more detail the statistics of the flow for a particular simulation. Fig. 11 shows the various ways to view histograms of the potential vorticity for simulation D from the forcing experiment. At the top of the figure the joint histogram for potential vorticity and streamfunction is shown for all zonal distances and times on the y = 0 (or approximately the ψ = 0 line). The histogram is trimodal, looking like a mountain range with three peaks running in a fairly straight q − ψ line. We see here three distinct statistical populations in q and ψ. Below in Fig. 11 are histograms of potential vorticity using various ways of taking the statistics. On the left we have shown the histogram of potential vorticity for all zonal distances and times on the ψ = 0 line. This is equivalent to the joint histogram above where we have marginalized out the streamfunction. We see that although the marginalized histogram has three peaks the joint histogram tells us that each peak corresponds to different values of instantaneous streamfunction.
In what follows, it will be shown that the two side peaks correspond to the regions of homogenized potential vorticity either side of the potential vorticity step. This gives rise to the skewness/kurtosis behaviour considered in Section 5 . In the middle plots, we consider how moving to a frame of reference travelling at the wave-speed of the most dominant Rossby wave may affect the statistics. In this way we hope to remove much of the meandering of the jet that is due to the dominant Rossby wave passing through the channel. This yields the Galilean transformed timeaveraged streamfunction shown in Fig. 11 ( bottom-middle ). Despite changing the relative size of the three peaks discussed previously, this method of taking the statistics is insufficient to remove the trimodality of the histogram. If however, we use the instantaneous streamfunction ψ = 0 as the spatial ordinate for taking the statistics we completely remove the two side peaks in the distribution giving rise to an unimodal histogram. Thus by explicitly moving with the meandering jet we have drastically altered, and simplified, the structure of the potential vorticity distribution. This also proves that the two side peaks, apparent when the statistics is viewed from an Eulerian perspective, are due to the large scale meandering of the mixing barrier and are distributions of the regions of homogenized potential vorticity either side. This leaves only the potential vorticity statistics within the strong potential vorticity gradient itself as can be seen in Fig. 11 . Fig. 12 shows the first four moments considered in previous sections now as a function of the normalized instantaneous streamfunction. The behaviour of the mean is qualitatively unchanged from that observed in Section 4 , while the standard deviation displays major qualitative differences. For the lower forcing strengths A-E the standard deviation is suppressed around ψ = 0 while for the two highest forcing strengths the central standard deviation is enhanced. The variation between simulations in the skewness and kurtosis is very large and it is hard to discern a consistent behaviour. However it clear that characteristic behaviour indicative of a mixing barrier is now lost and the kurtosis is on the whole leptokurtic and does not show the low (close to K = 1 ) kurtosis behaviour when the statistics was taken on contours of timeaveraged streamfunction. The change of sign of skewness is no longer coincident with low kurtosis values. Thus, we have shown that this behaviour disappears when the mixing barrier is removed either through changing the physics or by moving into a frame of reference moving with the meanders of the jet.
At this point it is worth considering how these statistics, given in Fig. 12 , correspond to the prediction of equilibrium statistical mechanics theory. An equilibrium, energy-potential enstrophy, theory ( Salmon et al., 1976;Jung et al., 2006 ) predicts that the mean will be linear with respect to streamfunction, the standard deviation will be constant, the skewness will have value, S = 0 , and the kurtosis will have value, K = 3 . We see that these four criteria are approximately met in the centre of the jet for the lowest strength of forcing. This suggests agreement with the results of the equilibrium theory when forcing is weak and the meandering of the mixing barrier has been taken into account. This chimes with the study of Jung et al. (2006) in which they moved into the frame of reference of their large scale Rossby wave to observe their statistics. We suggest that the reason this procedure does not work here is the absence of any significant scale separation in contrast to the experiment of Jung et al. (2006) . This could be due to various factors including: strength of dissipation; length scale of forcing; and strength of the β-term.
The joint histograms agree with the interpretation of separate statistical populations either side of zero streamfunction and is consistent with presence of mixing barrier. Here we considered one example of the joint q − ψ statistics but will consider the histograms for other simulations further in Section 7 . We have shown that performing the statistical analysis in a quasi-Lagrangian frame of reference removes the multimodal nature of the distribution simplifying the statistics but at the cost of a more complicated time-averaging. In addition, the simulation with lowest forcing strength gives an approximate agreement with Jung et al. (2006) in both the quasi-Lagrangian frame and the reversed forcing statistics. It could be argued that as the strength of forcing is reduced the statistics tends to the (energy-enstrophy) equilibrium statistical mechanics prediction, but only when the mixing barrier is removed.

Discussion
In this discussion we will draw together the results from the previous sections by considering: joint distributions of potential vorticity and streamfunction; the concept of mixed distributions and their importance to describing mixing barriers; implications for stochastic parameterizations; and the application to the South-ern Ocean. In this way we are able to highlight some key implications of the results of this study.

Joint probability distributions
To highlight the statistical nature of mixing barriers it is fruitful to consider the q − ψ joint histograms for all the forcing strengths in simulations A-G, given in Fig. 13 . As with the example given in Fig. 11 we have plotted the joint counts of potential vorticity and streamfunction for 100 × 100 bins, width set by the range of the data. As mentioned, the distributions are a mixture of bimodal and trimodal; A and G are bimodal; B and F have a weak central peak; C, D and E are trimodal. As described previously we see the clear separation of the statistical populations consistent with the picture of a meandering potential vorticity step, but further to this we see interesting features of the statistics.
Firstly the question is raised of what causes the existence of the third peak seen in B-F, centered at 0. From examining the details of the flow evolution it appears that this central peak coincides with the emergence of a spatially narrow region in the core of the jet where the potential vorticity gradient is weakened coming to zero (or even a weakly negative value). This narrow homogenized region appears intermittently within the otherwise strong potential vorticity gradient associated with the channel Rossby wave. This can be seen in the snapshot from simulation D, Fig. 3 . While difficult to observe in noisy snapshots of potential vorticity, this becomes apparent when the potential vorticity is plotted against streamfunction as the streamfunction's latitudinal gradient (or zonal velocity) has its largest magnitudes near the centre of the jet, see Fig. 11 , D. Secondly we notice that the areas of high bin count are elongated in the streamfunction direction, giving the same value of mean potential vorticity for many values of streamfunction: apparent in the data from simulations C-G. This quite clearly demonstrates the potential vorticity homogenization paradigm discussed in, for example, Esler (2008) . Indeed, the potential vorticity step with adjacent regions of homogenization is made very apparent when the flow is viewed as a potential vorticity/streamfunction joint distribution.

Mixed distributions
We have shown that the statistics of potential vorticity in the presence of a mixing barrier are described by mixed distributions, a sum of distributions of random variables. Mixed distributions are distributions corresponding to random variables created by a random choice of other random variables. In other words, they are a convolution of discrete and continuous random variables. Here, the discrete random variable represents the regions separated by the mixing barrier and the continuous random variable represents the statistics of the homogenized regions. The full Eulerian statistics is then given by randomly selecting which side of the mixing barrier we are standing through a Bernoulli trial. The simple models (17) and (19) are examples of mixed distributions. This opens up the possibility of using Gaussian Mixture Modelling methods to analyse the statistics of meandering mixing barriers. These methods have been used in atmospheric dynamics ( Hannachi, 2007 ) to consider weather regimes and also, more recently, in oceanography ( Maze et al., 2017 ) to classify hydrographic profiles.
An important consequence of the mixed nature of these distributions is that they cannot be constructed from a maximization of entropy given knowledge of an arbitrary number of moments. If we assume that we know that the first N moments of a probability, then maximization of information entropy gives a probability distribution of the form Mixed distributions such as (17) and (19) cannot be inferred in this way through only knowledge of the measured moments. The maximum entropy distribution, (20) , can be used to approximate multimodal distributions such as (17) and (19) but this requires N to be very large (example calculations suggest N needs to be as high as 100), considerably beyond what is available from data. Essentially this comes from the difficulty of representing a bimodal distribution by the exponential of a Taylor expansion: a very large number of terms would need to be kept to capture two well separated peaks. The maximum entropy method of approximating probability distribution functions was successfully used in Porta Mana and Zanna (2014) to fit the diagnosed distribution of a stochastic parameterization where the distributions considered were unimodal; this method of fitting probability distribution functions would be inappropriate for use in the case considered in this study.

Stochastic parameterization
The results in this study have important consequences for the way we approach the stochastic parameterization of ocean mesoscale turbulence. If we parameterize the Eulerian eddy forcing tendency, then it is essential that we do this in a way that will lead to a mixed probability distribution for potential vorticity in the presence of a strong mixing barrier. It is therefore prudent to both take into account higher order cumulants as well as spatial correlations as these are required to describe the structure of jets. This is the approach followed in stochastic structural stability theory (e.g. Farrell and Ioannou, 2003 ) and in the cumulant expansion method (e.g. Marston et al., 2008 ); by solving prognostic equations for statistical cumulants of the flow it is possible represent eddy-mean flow interactions which respect the steady state structure and dynamics of turbulent jets.
An alternative interpretation of the results presented in this study is that a stochastic parameterization should be implemented in a Lagrangian sense, allowing the noise to be simpler in structure: monomodal instead of the multimodal statistics seen in a Eulerian frame. Recently parameterizations have been developed which modify the Lagrangian tendencies rather than the Eulerian tendencies ( Holm and Nadiga, 2003;Holm and Wingate, 2005;Porta Mana and Zanna, 2014;Zanna et al., 2017 ); Porta Mana and Zanna (2014) do this by stochastically modifying the Lagrangian tendency of potential vorticity in a quasi-geostrophic model. We can see this by considering the potential vorticity equation

Dq Dt
where D, F and η are the dissipation, forcing and stochastic noise respectively. In this formulation the noise can only influence the Eulerian tendency of potential vorticity. However, the parameterization introduced by Porta Mana and Zanna (2014) can be written as Dq Dt where κ is now a stochastic noise with a negative mean value. By rewriting as we can clearly see that parameterizations of this form modify the Lagrangian tendency. This makes parameterizations of this form an advantageous framework for a stochastic treatment as the potential vorticity distributions become simpler when viewed in a Lagrangian frame of reference, see Section 6 .

Skewness and kurtosis in the Southern Ocean
So far, we have only considered the statistics of potential vorticity in a highly idealized barotropic model; now we consider the application of the ideas developed here to the Southern Ocean.
Due to the sparsity of in-situ observations, we calculate the statistics from the Southern Ocean State Estimate  ), a numerical model at 1/6 °resolution, constrained to fit, as closely as possible, a wide range of observational data. Hughes et al. (2010) showed that the skewness and kurtosis of relative vorticity can be used to identify strong jets from altimetry observations; using data from the Southern Ocean State Estimate we can look for this skewness/kurtosis signature at depth. To provide the best analogy to the barotropic model considered in this study as, ignoring diapycnal mixing, the potential vorticity is conserved on a neutral density surfaces with sources and sinks appearing only when the surface outcrops or interacts with topography. We will take statistics on 5-day mean neutral density surfaces ( Jackett and McDougall, 1997 ) of the Ertel potential vorticity where ω is the relative vorticity; f is the planetary vorticity; γ n is the neutral density; and ρ 0 is the reference density under the Boussinesq approximation. By calculating the potential vorticity from Southern Ocean State Estimate data and interpolating onto 5-day mean neutral density surfaces we are able to compute the skewness and kurtosis as functions of longitude, latitude and neutral density. Maps of skewness and kurtosis are given for two neutral density surfaces in Fig. 14 . Following the diagnostic tool, and intuition, developed in Section 5 we look for a positive-negative feature in the skewness concomitant with a high-low-high feature in the kurtosis. Fig. 14 a and b have several features of these type (highlighted in boxes) and correspond to regions of high potential vorticity gradient within the Antarctic Circumpolar Current as well as the Malvinas current (east of South America). Although weaker in other places the skewness/kurtosis feature seems to extend through the whole domain. Looking at Fig. 14 c and d, this is more apparent; now a connected feature can be seen circumnavigating the Southern Ocean. It is interesting to note that the feature in Fig. 14 c and d, has shifted south compared to Fig. 14 a and b; the southward migration of the skewness/kurtosis feature be seen very clearly when additional density surfaces are examined. Explaining this migration of the jets with increasing density lies beyond the scope of the present study but this migration underlines the usefulness of this framework in studying the vertical structure of jets in the Southern Ocean using primitive equation ocean models.
We have shown that the skewness/kurtosis signature provides a diagnostic tool with which mixing barriers can be correctly identified and will be important for future studies on the influence of these jets on Southern Ocean mixing.

Concluding remarks
The relationships between streamfunction and the various statistics we have diagnosed in this study cannot be explained by equilibrium energy-potential enstrophy statistical mechanics theories, contrary to the result of Jung et al. (2006) who found good agreement with a barotropic rotating annulus experiment. Whether more general versions of the theory can be used to explain the generic statistical characteristics of this flow will require a deeper understanding of the response of the statistics to forcing and dissipation. In this study we explored the phenomenology of these statistics over a wide range of forcing and dissipation strengths. We have shown that the mean and standard deviation of the potential vorticity as a function of time-averaged streamfunction exhibits self-similar behaviour for barotropic flow across a wide magnitude range of eastward wind stress. The self-similarity is also observed when the dissipation is perturbed and when the direction of wind stress is reversed, although in these cases selfsimilarity is less clear. As the self-similarity is stronger for eastward forcing, and consequently eastward mean flow, we infer that this self-similar behaviour is influenced by the emergent mixing barrier that forms in the present set of numerical experiments. Additionally, by relating the scaling factor of this self-similarity to the global potential enstrophy of the flow, we have related the dependence of the mean and standard deviation of potential vorticity to the global balanced quantities of the flow in a statistically steady state. Importantly, this reduces the problem of determining the mean and standard deviation to only needing to know a single quantity: the potential enstrophy (or eddy potential enstrophy in the case of the standard deviation). Further work will be needed to fully determine the range of parameter space for which this self similarity behaviour exists. It is also clear that the implications of this result have not been fully elucidated with respect to eddy-mean flow interaction in a statistically steady state jet. In future work it is important to relate this self-similar scaling to the imposed parameters of the model such as the wind stress strength.
The behaviour of skewness and kurtosis as a function of timeaveraged streamfunction is characteristic of a mixing barrier and a simple statistical model was used to paint a statistical picture of a mixing barrier. Any potential stochastic parameterization must respect the multimodal, and mixed, statistics of a mixing barrier along with its consequent skewness and kurtosis behaviour. This highlights the importance of understanding the influence of highly non-Gaussian noise on the mean flow. It has been shown that the statistics is no longer multimodal when viewed from a frame of reference moving with the barrier itself. This suggests a separation between the statistics due to the diffusive eddies in the homogenized regions of potential vorticity and the statistics of the barrier meandering, that is, the statistics of the large scale Rossby waves. This also highlights the potential importance of Lagrangian thinking when attempting to parameterize eddy statistics.
We have argued that the results presented here have important implications for stochastic parameterizations as well as being an useful tool for the study of ocean jets and their influence on mixing. Despite the increased complexity of comprehensive ocean models, the statistical picture we have painted relies on the existence of a mixing barrier rather than on the barotropic nature of the dynamics. For this reason, we anticipate that our results should have much broader applicability to baroclinic flows (cf. Hughes et al., 2010 ), and seen here using Southern Ocean State Estimate data. In future development of this work, it will be important to relate the neutral density dependence of the statistics to the depth dependence of the statistics; we expect strong mixing barriers forming near the surface but not at depth, consistent with observations of the Gulf Stream ( Bower et al., 1985 ), observationconstrained data in the Southern Ocean ( Abernathey et al., 2010 ), and idealized numerical calculations ( Esler, 2008 ).
In conclusion, we have presented a simple yet useful description of the statistics of a barotropic jet with a mixing barrier over a wide range of forcing and dissipation strengths. The jet meandering can be described by Bernoulli statistics while the diffusive turbulence either side is described by a superposed small-scale noise. This way in which the mixing barrier separates the statistical populations of potential vorticity either side can be considered not only as a statistical model of mixing barriers but as a expression of what a mixing barrier is in a statistical sense. We hope that the statistical picture of mixing barriers presented here can provide a foundation for future work on the stochastic representation of eddies and their dynamics in ocean models.