Properties of the Subantarctic Front and Polar Front from the skewness of sea level anomaly

The region of the Southern Ocean that encompasses the Subantarctic Front (SAF) to the north and the Polar Front (PF) to the south contains most of the transport of the Antarctic Circumpolar Current. Here skewness of sea level anomaly (SLA) from 1992 to 2013 is coupled with a meandering Gaussian jetw model to estimate the mean position, meridional width, and the percent variance that each front contrib-utes to total SLA variability. The SAF and PF have comparable widths (85 km) in the circumpolar average, but their widths differ signiﬁcantly in the East Paciﬁc Basin (85 and 60 km, respectively). Interannual variability in the positions of the SAF and PF are also estimated using annual subsets of the SLA data from 1993 to 2012. The PF position has enhanced variability near strong topographic features such as the Kerguelen Plateau, the Campbell Plateau east of New Zealand, and downstream of Drake Passage. Neither the SAF nor the PF showed a robust meridional trend over the 20 year period. The Southern Annular Mode was signiﬁcantly correlated with basin-averaged SAF and PF positions in the East Paciﬁc and with the PF south of Australia. A correlation between the PF and the basin-scale wind stress curl anomaly was also found in the western extratropical Paciﬁc but not in other basins.


Introduction
The Antarctic Circumpolar Current (ACC) is the one of the most dominant, dynamical features in the Southern Ocean. It connects the Southern Hemisphere Pacific, Atlantic, and Indian Ocean basins. Much of the large-scale, eastward flow organizes into jets in three distinct regions (from north to south): the Subantarctic Front (SAF), the Polar Front (PF), and the Southern ACC Front (SACCF). Each of these fronts characterizes regions where near-surface hydrographic properties change rapidly. The fronts tend to be aligned zonally except in regions where steeply sloping topography steers their paths, for example south of New Zealand around the Campbell Plateau, Scotia Sea and in the Brazil Basin. To first order, these fronts are in geostrophic balance resulting in jets whose zonal velocities are described by u g 52 g f @g @y ; (1) where u g is the zonal geostrophic velocity, g is gravitational acceleration, f is the Coriolis frequency, and g is the geopotential height and is equivalent to the sum of the mean dynamic topography (MDT) and sea level anomaly (SLA) at any given time. The circumpolar path of these fronts was estimated using hydrographic data to define criteria describing each of these fronts yielding the traditional view of three coherent, circumpolar structures [Orsi et al., 1995]. Subsequent studies found evidence that the SAF and PF split into multiple branches in some regions, suggesting a more complicated frontal structure [Belkin and Gordon, 1996;Sokolov and Rintoul, 2002].
ensure that any given maximum is related to a jet. Additionally, these methods also require an a priori choice of a threshold value.
A commonly used method of determining the path of ACC fronts relies on estimates of MDT and is based on the interpretation of SSH contours as geostrophic streamlines. One of the earliest applications this ''contour method'' matched density fronts from hydrographic and expendable bathythermographs to SSH contours from satellite altimetry south of Tasmania [Sokolov and Rintoul, 2002]. They showed that these multiple branches sometimes merge to form a single strong jet or diverge to yield two weaker, but distinct, cores. Marrying the concepts of hydrographic fronts with dynamic jets, Sokolov and Rintoul [2002] further noted that the characteristic SSH contours of these fronts were coincident with large SSH gradients and enhanced geostrophic velocities via equation (1). Sokolov and Rintoul [2007] characterized fronts as being the SSH contours that best fit high gradients of SSH when integrated over its entire circumpolar path. Under the assumption that these SSH contours represent the same front over time, this method can be used to track the path of these fronts over time [Sall ee et al., 2008;Rintoul, 2007, 2009a; Kim and Orsi, 2014].
Recent studies have evaluated this assumption, questioning whether the SSH contours associated with the SAF and PF are actually continuously circumpolar fronts [Graham et al., 2012;Chapman, 2014] and also asking whether these contours could shift in response to steric changes in the large-scale sea level field [Gille, 2014]. Graham et al. [2012] and Chapman [2014] note too that the process of ascribing a contour of SSH by maximizing the circumpolar integration of SSH gradient does not guarantee that any location along the path is a local maximum in SSH gradient. Finally, the contour method yields different frontal locations and number of fronts depending on the SSH data product used as they can differ in both their estimates of mean SSH and SSH gradient [Volkov and Zlotnicki, 2012].
The method detailed in this study is based on the findings of Thompson and Demirov [2006;hereafter TD06] that meandering jet-like features are associated with non-Gaussian probability density functions (PDFs) of SLA. The shape of this distribution can be well characterized by skewness and kurtosis [Hughes et al., 2010] and have been used to separate jet-like features from eddies [Chapman, 2014]. Skewness and kurtosis are defined as the third and fourth moments about the mean, normalized by the standard deviation: where by convention n 5 1 for skewness and n 5 2 for kurtosis, X is the distribution of SLA at an individual location, E is the expectation operator for a distribution of quantities v with N elements v i defined as X N i50 v i =N, l is the mean of these SLA, and r is the standard deviation. TD06 show in a quasi-geostrophic model that skewness is positive poleward of a jet and negative equatorward. Hughes et al. [2010] find that, when used together, skewness and kurtosis describe most of the spatial variability in PDFs of SLA in the ocean They further demonstrate that jet cores are associated with PDFs with zero skewness and low kurtosis.
In section 2, we extend TD06's findings using a meandering Gaussian jet model to better understand how different properties of a jet can modify the magnitude of skewness. In section 3, we apply this method to the Southern Ocean to identify the mean position, width and local effect on SLA variability of the SAF and PF. In section 4, we subset the data annually to determine whether a trend exists in the meridional position of fronts, and to investigate whether the interannual variability in the frontal position can be ascribed to local variability in forcing or to two large patterns of climate variability: the El-Niño Southern Oscillation (ENSO) and the Southern Annular Mode (SAM). Section 5 discusses our findings in the context of studies of SAF and PF variability and compares the benefits of our front detection method relative to others. Finally, section 6 summarizes the main findings of this study.

Skewness of SLA
The skewness statistic is a measure of the symmetry of a distribution whose values are largely determined by the tails of the distribution. For example, a perfectly symmetric distribution (e.g., a Gaussian) has zero skewness, a distribution whose right tail is thicker than its left tail (e.g., an Inverse Gaussian) has positive Journal of Geophysical Research: Oceans 10.1002/2015JC010723 skewness, and a distribution whose left tail is thicker than its right tail has negative skewness. TD06 showed that the presences of eddies as well as strong geostrophic currents can cause strong nonzero regions of skewness. An eddy field dominated by either cyclonic or anticyclonic rotation can result in an SLA distribution with negative or positive skewness, respectively. Using an idealized quasi-geostrophic model and skewness calculated from satellite altimetry, TD06 further demonstrated that oceanic jets have a distinct signature where skewness transitions from positive south of the jet to negative north of a jet. South of a jet, SSH is low, whereas SSH is high north of a jet. If a jet meanders southward, some anomalously high values of SSH will be included in the distribution resulting in a PDF with positive skewness (Figure 1d). An analogous argument explains why skewness is positive poleward of a jet (Figure 1b).

The Gaussian Meandering Jet Model 2.2.1. Description of the Idealized Jet Model
To be able to investigate the location and size of the skewness transition, we consider an idealized model for a meandering jet (Figure 1) which has a purely zonal Gaussian velocity profile and meanders in time according to u g ðy; tÞ5u 0 exp 2 y2 y 1dðtÞ where y is the distance from the jet center, t is a time variable, u 0 is the maximal zonal velocity (assumed to be constant in time), y is the mean position of the jet, d(t) is a normal distribution with zero mean and a°°F igure 1. Synthetic skewness curves generated using a meandering jet centered at 608 S. (a) Ten realizations of SSH using the Monte Carlo method described in section 2, (b-d) the probability density function of sea surface height south of the jet (618 S), at the center of the jet (608 S), and north of the jet (598 S), and (e) the skewness of SLA at each latitude demonstrating the transition from positive to negative skewness as the jet is traversed from south to north.
Journal of Geophysical Research: Oceans 10.1002/2015JC010723 standard deviation r which represents the distance over which the jet meanders, and L is the width of the jet [Gille, 1994]. Using the geostrophic relationship (equation (1)), the corresponding SSH will be g g ðy; tÞ5g 0 erf y2 y 1dðtÞ ffiffi ffi 2 p L ; where g 0 is equal to p=2 ð Þ 1=2 fLu 0 =g and can be interpreted as the size of the SSH jump associated with the jet. Averaging in time yields the equivalent of mean dynamic topography, g g ðyÞ5g 0 erf ðy2 yÞ= ffiffi ffi 2 p L À Á . Here g g ðy; tÞ represents an ocean in which all variability is associated with meandering currents, whereas the Southern Ocean also contains a large number of coherent eddies [e.g., Keffer and Holloway, 1988;Fu, 2009]. The total SLA (e.g., as observed by a satellite) from both the meandering jet and the eddy field is then g obs ðy; tÞ5g g ðy; tÞ1A / /ðy; tÞ; where A / is the amplitude of the sea level anomaly from the eddy field and /ðy; tÞ is a structure function representing the SLA of the eddy field, which we assume to be approximated by white noise (discussed in further detail in the next section).
Despite / having a mean of zero, the relative contribution of g g and A / to total SLA (g obs ) does affect the skewness, because if A / is larger than g 0 the SLA distribution is dominated by noise rather than the SLA of the jet. We now define a new term SNR5g 0 =g obs that represents the fractional contribution of the jet signal to the total observed SLA versus the eddy ''noise'' A / . If SNR is equal to one, the observed SLA is due solely to the meandering jet in contrast to the scenario where the total SLA is dominated by eddy ''noise'' (SNR 5 0). Here both A / and g 0 are assumed to be constant in time such that the maximum and minimum g obs at any point in time is similar. Additionally, as can be shown from equation (2), for a distribution X, c 1 ðXÞ5c 1 ðaXÞ where a is a scalar. This allows equation (5) to be rewritten using the SNR term to express SLA as the sum of the fractional contribution of the jet and the eddy field without changing the magnitude of skewness g N ðy; tÞ5SNR erf y2 y1dðtÞ 2 ffiffi ffi 2 p L ! 1ð12SNRÞ/ðy; tÞ; where g N ranges from 21 to 1 if SNR 5 1. Using results from Gille [1994], the meander distance r is set to 1.7 times the width of the jet. Each of the model parameters have been summarized in Table 1 for reference.
A skewness curve can be generated using a Monte-Carlo simulation by choosing values for the three free parameters y, L, and SNR and drawing 1000 samples from d(t). This large number of samples helps ensure that the random numbers used as an analog for the meandering of the jet and for eddy variability have zero skewness. An example for a purely zonal, meandering, Gaussian jet in the Southern Ocean with y560 S, L 5 80 km, r 5 136 km, and SNR 5 0.8 is shown in Figure 1 Figure 1d) are skewed while the PDF at the center of the jet, 608S has close to zero skewness ( Figure 1c).

Relationship Between Skewness and Model Parameters
Plots of SLA skewness in TD06 were normalized such that the maximum or minimum of the skewness transition was equal to 1. In the example (Figure 1), skewness reaches a maximum of around 1.5. Here we vary the shape parameters in equation (6) individually to gain insight into how the properties of the jet can affect the skewness of the resulting SLA distribution. The mean position parameter y has no effect on the magnitude of the skewness, but does shift the zero crossing of the skewness curve (not shown). The effects of the other two shape parameters, L and SNR, on the skewness transition were determined empirically. To investigate L, the same parameters as those shown in Figure 1 were used but L was allowed to vary from 20 to 120 km. Similarly to investigate SNR, L was set fixed to 80 km, but SNR was allowed to vary from 0.01 to 1. This empirical sensitivity test shows that the maximum skewness is inversely proportional to the jet width L ( Figure 2a) and that SNR has the largest impact on the maximum skewness. An SNR of approximately 0.6 is required for the signal of the jet to be distinguishable from eddy variability ( Figure 2b). Additionally, because SNR and maximum skewness are so closely related, strong jets which dominate SLA variability are readily identified as regions with large skewness transitions.
In the examples above and for the remainder of this study, we assume that the eddy SLA signature (/ in equation (6)) can be modeled as white noise. This is potentially problematic as eddy variability has been estimated to be spatially correlated on length scales up to about 85 km [Gille and Kelly, 1996]. To test the effect of this assumption, we performed the same sensitivity analysis using colored noise with spectral slopes from 22 to 25 (which covers the range reported in the literature) [e.g., Xu and Fu, 2011;McCaffrey, 2014] and a wavelength roll-off beginning at 85 km. Both the magnitude of skewness ( Figure 2) and the shape of the skewness curve (not shown) were affected minimally by the spectral content of the noise. In the case where L varied with a fixed SNR, the maximum skewness with white noise was 5% lower than with colored noise. Sampling multiple realizations of the colored noise in time yielded Gaussian PDFs, implying that the spatially distributed colored noise when sampled in time behaves like white noise.
Jets may be in close enough proximity to each other such that their velocity envelopes overlap [Gille, 1994].
To investigate the impact on skewness of two jets in close proximity, we simulate this case under the assumption that the jet-jet interactions are linear (not shown). Skewness is a nonlinear statistic and so the individual skewness profile of each jet cannot be linearly summed [Hughes et al., 2010, Appendix A]. Here we find that for two jets a and b, if the distance between their centers j y a 0 2 y b 0 j is less than the sum of their widths L a 1 L b , only one skewness transition will be observed, asymmetrically centered around the peak of the summed velocity. The zero crossing of this skewness transition still corresponds roughly to the average of the two jet centers. When applied to observed skewness, this implies that given an asymmetric skewness transition the zero crossing is still a robust estimate of jet locations. However, in regions where multiple jets overlap, our use of a single Gaussian may yield overly broad estimates of the jet width.

Application to the Southern Ocean
To determine ACC jet positions from observations, we use the delayed-time, reference along-track SLA product provided by AVISO. This product represents a continuous measurement record from the TOPEX/ Poseidon, Jason-1, and Jason-2 satellite missions. Each satellite measures along a fixed reference orbit with . The effect of changing the shape parameters of the meandering Gaussian jet model on the magnitude of the resulting skewness using red noise (blue) and white noise (black). (a) L ranges from 20 to 120 km (SNR 5 0.8, y 5 0) and (b) SNR ranges from 0 to 1 (L 5 80 km, y 50). See Table 1 and text for explanation of parameters.
Journal of Geophysical Research: Oceans 10.1002/2015JC010723 254 ground tracks (127 descending and 127 ascending). SLA measurements are repeated every 9.9 days from 25 September 1992 to 7 August 2013 with an approximate 10 km along-track resolution. Unfiltered along-track data are used for this study rather than a gridded product, because the gridding process affects higher moment statistics such as skewness [Sura and Gille, 2010]. Additionally, in the Southern Ocean the first baroclinic Rossby radius of deformation ranges between 10 and 30 km [Chelton et al., 2011] and is better resolved by the along-track data than by the gridded product. The added measurement noise associated with the unfiltered data is assumed to be white noise and is incorporated as part of the A / term. Following TD06, we calculate skewness of SLA along every track to show regions in the Southern Ocean where skewness transitions from positive to negative. Throughout this study, skewness was calculated using a biased estimator (differences using an unbiased estimator were negligible). The positions of the middle branches of the SAF and PF estimated by Sokolov and Rintoul [2009b] (red and blue lines, Figure 3) roughly follow regions where skewness transitions between positive (red) and negative (blue) values. The along-track skewness of SLA from the entire altimeter record on ground track 6 (location shown in Figure 3) shows these distinct skewness transitions more clearly (the two largest indicated by the blue and red colored portions of the skewness curve in Figure 4). The skewness is fairly noisy on this track (and on all others) which leads to false positives when trying to identify skewness transitions by eye. In section 4, the skewness of SLA was calculated based on annual subsets of the data between 1993 and 2012, when measurements of SLA were collected over the entire year. To evaluate the sensitivity of our results to using calendar year subsets of the data, we perform the same analyses over year beginning in April and ending in March (which is more representative of ENSO) [Stein et al., 2010], which will be referred to as the ENSO-year subset. If frontal positions exhibit any statistically significant trends or  correlations to climate modes of variability in both subsets of the data, the result is more likely to be robust and nonspurious.
The following algorithm is applied to each track in order to consistently identify skewness transitions associated with either the SAF or the PF. First, the latitude bounds in which we look for skewness transitions are chosen to restrict the analysis to the ACC. Here we use the Sokolov and Rintoul [2009b] estimates of the middle branches of the SAF and PF as the basis for the envelope. Their estimates of front location are also based on satellite altimetry and provide higher spatial resolution than the hydrographic estimates of Orsi et al. [1995], though our results should not be sensitive to this choice, as the estimates from Sokolov and Rintoul [2009b] broadly agree with those of Orsi et al. [1995]. The northern bound of our search envelope is 58 equatorward of the SAF-M, and the southern bound is 58 poleward of the PF-M. Second, either the maximum or the minimum of a skewness transition must be statistically significant. In this study, statistical significance is determined using a bootstrap method [Efron and Tibshirani, 1994] at the 95% confidence level. Alternatively, a significance level can be calculated directly using the standard error of skewness: , where n is the number of samples [Fisher, 1930]. Finally, because the largest skewness transitions should be associated with strong jets (i.e., higher SNR), the two largest skewness transitions along each track are chosen. In this study, the northern transition is referred to as the SAF and the southern one is referred to as the PF. Note here that these labels are used because these transitions tend to roughly align with the SAF and PF locations from Sokolov and Rintoul [2009b]. However, in the absence of any other information (e.g., SSH, SST, or subsurface hydrographic gradients), these frontal locations may actually correspond to other strong fronts such as the SACCF, the Subtropical Front [Graham and De Boer, 2013], or a different branch of the same front.
For each of these skewness transitions, the shape parameters y, L, and SNR are determined by a nonlinear least squares fitting routine minimizing the normalized sum of squared residuals (NSSR) ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðc obs 2c s Þ 2 =ðN23Þ q , where c obs is the observed skewness, c s is the skewness derived from the Monte-Carlo simulation in section 2.2, N is the number of data points used to fit, and the 3 arises from the number of free model parameters. Based on examining a number of fitted curves, we define a good fit where NSSR is < 0.2. This fitting method relies on randomly drawn samples and so is susceptible to uncertainty in the shape parameters. To estimate the inherent error in this method, we generated 100 independent skewness curves using the Monte-Carlo method described in section 2 with L 5 100 km, y5250 , and SNR 5 0.8 and evaluated if the fitting routine can reproduce the parameters. When the number of simulated measurements is large (here 700, a number similar to the full data record), the difference between the actual and estimated values for L was 20 6 1 km, y was 0:160:01 , and SNR was 0.08 6 0.01. However, when applied to smaller subsets (i.e., 36 measurements, similar to the size of the annual distributions of SLA), the mean position of the front, y, is still well constrained, but L and SNR are not. This behavior arises from the sensitivity of skewness to outliers. For a large number of samples, this effect is suppressed, but for smaller realizations, the Monte-Carlo method used here does not consistently generate the same width or magnitude of the skewness transition. Finally, the width estimates of the jet may be biased because satellite tracks are not always perpendicular to the axis of the jet. In this case, the width that a satellite ''sees'' may be artificially broadened, instead measuring the along-track width.

Long-Term Mean Position and Width
As an initial investigation of the accuracy of the method, we compare the estimated mean jet positions from our method with those determined by the frontal positions from Sokolov and Rintoul [2009b]. In general, consistent with previous studies [e.g., Orsi et al., 1995;Graham et al., 2012;Sokolov and Rintoul, 2009b], we find that frontal locations align along zonal ridges and are steered around areas of significant topographic features (Figure 5a). Dynamically, this topographic control is due to near conservation of potential vorticity (PV) along the path of a fluid parcel which in the barotropic framework is f/H, where f is the Coriolis frequency, and H is the height of the fluid [Marshall, 1995]. If the fluid encounters a ridge, H decreases, so f must also decrease to conserve PV, implying an equatorward deflection.
South of Africa between 08E and 308E, the SAF (red dots in Figure 5a) flows across a small abyssal plain before reaching the Southwestern Indian Ridge (308E-508E) where the front is primarily constrained along The PF (blue dots in Figure 5a) follows topographic features similar to those followed by the SAF. At the Southwest Indian Ridge, the PF and SAF are located within a few degrees of each other. The PF then sharply diverges from the SAF at Kerguelen Plateau southward through the Fawn Trough south of the Kerguelen Islands. It then tracks the crest of the Southeast Indian Ridge closely before traversing it at around 2208E and flows south along the southern part of the Bellingshausen Basin. The PF and SAF begin to converge through Drake Passage (but remain distinct from each other) until they traverse the Mid-Atlantic Ridge (around 3408E).
The width of both the SAF (87 6 10 km) and PF (80 6 8 km) are similar when averaged over their circumpolar paths. As a caveat when comparing these to other estimates of jet widths in the ACC, these quantities refer to the along track, Gaussian width of each jet's velocity profile. Our estimates are almost twice as wide those in Gille [1994] potentially because of our assumption that each front represents a single Gaussian jet. As discussed in section 2.2.2, two jets whose velocity envelopes overlap appear as a single, broad skewness transition.
The largest differences in width between the SAF and PF occur when their paths are the furthest apart (1808E and 2408E-3008E) where the PF narrows to around 60 km. The narrowing of the PF may be a result of changes in the baroclinic Rossby radius [Chelton et al., 2011] which decreases poleward. In addition, topography appears to influence the width of the SAF and PF ( Figure 5). In the region of the Kerguelen Plateau (608E-808E), the PF width decreases from 90 to about 75 km and the SAF narrows downstream (808E)°°°°°°°°°°°°°F  Figure 5a are the same shown in Figure 3. Topography shallower than 2000 m [Smith, 1997] is denoted by the gray shaded region.

10.1002/2015JC010723
from around 100 to 80 km. The SAF also narrows after passing the Southeast Indian Ridge. When both fronts reach the constriction between the Campbell Plateau and the Southeast Indian Ridge (1708E), the flow narrows considerably with the width of the PF decreasing from 90 to nearly 50 km and the SAF decreasing from 100 to 70 km. The SAF broadens over the Bellingshausen Basin after the Pacific Antarctic Ridge, whereas the PF reaches its narrowest width (other than the Campbell Plateau). The larger spatial variability of the PF's width and position compared to the SAF may be related to the PF's proximity to strong topographic features. For example, in regions near strong topography, the branches of the PF are restricted to a smaller meridional region which may sharpen the jets themselves or reduce the distance between jet cores, both of which would result in narrower widths.
The SNR of the PF anticorrelates with its width, suggesting that when the jet narrows because of topographic meridional constraints, the PF tends to have a stronger signature in SLA. In these regions, topography tends to constrain the flow to organize into a smaller number of jets. If the overall transport does not change in this region, then the velocities of these jets must increase, which would yield a higher signature in SLA [Graham et al., 2012]. Similarly, the SNR of the SAF decreases over abyssal plains (e.g., 1808-2008 and 2408-3008) where the lack of topography permits the SAF to become a more diffuse feature or to split into multiple jets [Thompson and Sall ee, 2012]. In the Indian Ocean (308E-1208E), both the SAF and the PF tend to have relatively low SNRs despite the strong topography in that basin which may be due to eddies spawned from the Agulhas Current and enhanced baroclinic instability downstream of the Kerguelen Plateau .

Temporal Variability in Frontal Locations
The long-term features of the SAF and PF (i.e., the mean path, width, and SNR) seem to be heavily controlled by bathymetry. However, on shorter time scales, the ACC may have significant variability caused by changes in the Southern Ocean westerly winds [Allison et al., 2010] or that may be internally generated variability in the ocean [Hogg and Blundell, 2006]. To determine how much the SAF and PF vary interannually, the mean latitude of each front is calculated by averaging the centers determined from the annual subsets of SLA over 58 longitudinal zones (centered at 2.58E, 7.58E, etc.).
As a simple measure of the variability in the locations of these two fronts, the standard deviation (r y ) of the annual frontal position in each longitudinal sector between 1993 and 2012 was calculated. The SAF's r y (red line in Figure 6a) remains fairly constant geographically with an average of 110 6 40 km. In contrast, the PF (blue line in Figure 6a) generally has much higher variability (170 6 57 km) and shows four distinct regions where r y is large: (i) the Brazil Basin (3008E), (ii) in the Bellingshausen Basin (2608E), (iii) along the Southeast Indian Ridge (1908E), (iv) upstream of Kerguelen Plateau at Elan Bank (608E). The elevated peak in region (i) is relatively narrow and may be an artifact of our methodology, which may be picking a more southern branch of the SAF or the PF in some years. In region (ii), the PF flows over an abyssal plain which leaves the jet free to wander over a wide meridional range without modifying its PV. In regions (iii) and (iv), the PF either flows very near to or over topographic obstruction, suggesting a relationship between variability in jet positions and topography, an interaction discussed by Thompson et al. [2010] and Thompson and Richards [2011]. Strong topography causes large changes in PV, and interactions with the mean-flow can cause a jet to drift toward low PV gradient regions [see Thompson and Richards, 2011, section 2.1 for more details] before the jet eventually dissipates. The peak in PF variability at 1908E might be explained by this mechanism, with the PF being close to the ridge in some years and further northward in other years as it migrates away from the sharp topographic gradient. At (iv), Elan Bank represents a bathymetric feature around which the flow may bifurcate. The enhanced variability here may then be related to the ''jet jumping'' mechanism described by Chapman and Hogg [2013] and Chapman and Morrow [2014] in which transport can transition from one jet to another on yearly or slower time scales. As our method would preferentially choose the stronger jet, this would result in higher variability in the mean position of the PF.
In addition to the frontal meridional variability, two other potential features of interannual variability in the ACC are considered: long-term meridional trends and changes in meridional position in response to climate modes of variability. We evaluate the significance of trends and correlations at a 90% significance level of (p < 0.1).

10.1002/2015JC010723
The poleward trend in wind stress related to the contraction of the polar vortex [Marshall, 2003] has led to the hypothesis that the ACC might respond with a similar poleward shift. Evidence for this exists in some observational studies [Sall ee et al., 2008;Sokolov and Rintoul, 2009a] and some low-resolution modeling studies [Saenko et al., 2005;Fyfe and Saenko, 2006] but there are a number of studies that suggest that the ACC is relatively insensitive to changes in wind stress [e.g., Gille, 2014;Meijers et al., 2012;Graham et al., 2012] (see also section 5 for further details). Here we find that neither the SAF nor the PF has a statistically significant trend averaged over their entire circumpolar paths. On smaller spatial scales, the PF has no statistically significant trends (blue dots in Figure 6b) whereas the SAF displays some zones with statistically significant trends (red stems in Figure 6b). Just eight of the 58 wide zones have statistically significant trends. These are likely spurious false positives considering that zones neighboring the significant regions are often of the opposite sign, and the trends lack any coherent, zonal structure. Additionally, none of these same zones show statistically significant trends in the ENSO-year subset of the data. These results suggest that a robust signature of a trend cannot be found in any basin.
SAM is the dominant mode of atmospheric variability in the Southern Hemisphere and can modify surface wind stress, particularly in the eastern Indian Ocean, southwest Pacific Ocean, and (to a lesser extent) Atlantic Ocean [Gong and Wang, 1999]. ENSO is one of the leading modes of global climate variability with atmospheric teleconnections to the Southern Hemisphere that project strongly onto SAM [Cai et al., 2011;Ciasto and Thompson, 2008]. Indices of the strength of SAM and ENSO are taken from those made available by the Ocean Observations Panel for Climate (http://ioc-goos-oopc.org/). These indices were annually averaged to match the same annual time scales for which skewness was calculated.
Correlations of SAM with mean SAF and PF latitude in the 58 longitudinal zones yield only a handful of statistically significant locations ( Figure 7a). As with the long-term trends, the significant correlations between SAM and the SAF position tend not to be spatially coherent except potentially in the southeast Indian Ocean (1208-1808). This may be due to the combination of SAM being particularly strong in this region and a relative lack of topography to the north of the jet allowing the SAF to shift in response to large SAM events. However, once again these correlations might be false positives as no two neighboring zones show significant correlations to SAM. For example, the zones centered at 142.58E and at 152.58E are of opposite sign. No coherent signal can be seen between SAM and the PF. Likewise, the statistically significant correlations between ENSO and these fronts are too random to indicate a strong relationship (Figure 7b). Similar to the trend analysis, these zones are also not significant in the ENSOyear analysis.
To examine basin-scale variability in SAF and PF location, we divide the ocean into Atlantic (3008E-258E), West Indian (258E-1008E), south of Australia (1008E-1808E), West Pacific (1808 E-2308E), and East Pacific (2308E-3008E) zones. An index of variability for position of the SAF and PF in each ocean basin is calculated as the mean latitude of the frontal positions within each of those defined sectors. To examine whether local variability in forcing can cause variability in each basin, the first empirical orthogonal function (EOF) of monthly wind-stress curl (WSC) anomaly from the ERA-interim reanalysis [Dee et al., 2011] is calculated in each basin within 108 of the ACC. The first mode of this EOF is used as an index of variability. Every basin's EOF has a SAM-like response with a shift of the WSC southward for positive phases of the local index of variability.
These basin-scale averages of the position of the SAF and PF lead to higher correlations (Figure 8) compared to the track-by-track investigation. In the East Pacific, both the SAF and PF mean latitude correlate strongly (p < 0.05) with SAM, suggesting that both fronts shift slightly northward (the SAF by 39 km and the PF by 20 km during 1r positive SAM year). This is also the region where topographic features may have little influence because of a relatively large abyssal plain. SAM is also anticorrelated with the PF south of Australia, where the influence of SAM is particularly strong, representing a shift southward of the PF of about 41 km for a 1r SAM event. ENSO does not correlate significantly with either the PF or SAF in any basin, though the strongest correlation occurs with the SAF in the East Pacific where ENSO originates. In the West Pacific, wind-stress curl is the only significantly correlated forcing of the PF. Despite the strong influence of SAM here, the first mode EOF of WSC is not significantly correlated to the SAM index. The lack of correlations in most basins suggests that the relationship between the SAF and PF latitudinal position may not be directly related to large-scale windforcing.

Discussion
The lack of long-term trends and the weak correlations with SAM and ENSO reported above is counter to the other studies [Sokolov and Rintoul, 2009a;Sall ee et al., 2008;Kim and Orsi, 2014] that use the SSH contour method to infer frontal position variability. Part of the disparity may arise from the fact that the results discussed here rely on annual estimates of frontal positions and ignore the possibility that the relationships to SAM occur on shorter or seasonal time scales. However, recent studies have noted a significant downside in using the contour method to infer temporal variability, notably that SSH contours may not be colocated with strong gradients in SSH along the entirety of their circumpolar path [Chapman, 2014;Graham et al., 2012]. In particular, Graham et al. [2012] argue that in these regions with weak gradients, SSH contours can migrate in time without necessarily corresponding to a change in transport which may lead to erroneous estimates of temporal variability.
Our finding that the path of the ACC is insensitive to changes in the large-scale wind field aligns with a number of other studies using independent methodologies. Gille [2014] diagnosed the meridional variability of the ACC by calculating the zonally averaged geostrophic transport-weighted latitude and found no long-term trend and only a weak correlation with SAM. In an analysis of historical model runs submitted to the fifth Coupled Model Intercomparison Project, the ACC was found to not shift significantly in response to a poleward migration of Southern Ocean westerlies [Meijers et al., 2012]. In a high-resolution coupled climate model configured for a 21st century warming scenario, Graham et al. [2012], showed that the position of the ACC remained virtually unchanged despite an overall shift in the westerlies of 1.38 poleward. Domingues et al. [2014] using a combination of in situ and satellite observations south of Africa found that the positions of the SAF and PF fronts were uncorrelated with SAM.
The methodology introduced in this paper is an alternative method to detecting fronts in the Southern Ocean (and perhaps elsewhere). It is relatively simple to implement, as skewness can be calculated directly from SLA without incorporating estimates of MDT. Additionally, unlike many other jet finding schemes that require an a priori estimate of the number of fronts to find or thresholding criteria [Chapman, 2014], this method can be modified to be entirely objective. Using the standard error of skewness, a significance threshold can be calculated and compared to the along track skewness (see Figure 4) to find statistically significant transitions. If only the mean position is needed, we found that extrapolating the zero skewness crossing using a linear fit to the skewness transition is a good approximation to the more complicated nonlinear fitting routine employed in this study (not shown, but note the roughly linear transition from positive to negative skewness in Figure 1a between 2628 and 2588). Also, because higher SNR corresponds with higher maximum skewness, jets can be ordered in terms of relative strength based on the magnitude of the transition.
This method has a number of limitations that could be improved upon in future studies. Other detection methods provide a picture of the ACC as being composed of closely spaced jets [Chapman, 2014;Sokolov and Rintoul, 2009a]. This suggests that perhaps the SAF or PF could be better approximated by a double Gaussian jet model [Gille, 1994;Hughes et al., 2010]; further investigation, however, must be done to determine whether this increases the number of free parameters beyond a tractable limit. Finally, particularly for studies that attempt to find every possible jet along a track, a kurtosis criterion could also be enforced when choosing skewness transitions, following the finding of [Hughes et al., 2010] that low kurtosis represents a barrier to mixing.

Conclusion
In this study, we used a Monte-Carlo simulation based on a meandering Gaussian jet model to demonstrate that a jet's width and SNR modulate the skewness signature first explored by TD06. Coupling this with a fitting routine allows for the robust identification of jets from the along track SLA data set. Unlike many other frontal detection schemes, no estimates of MDT are required, ensuring that any temporal variability in frontal location is unaffected by large-scale SSH changes. Also, this method requires no threshold criterion to be defined and any skewness transition that exceeds a statistical threshold potentially corresponds to a jet-like feature.
This method was then applied to the SLA in the Southern Ocean where the two largest skewness transitions along each satellite track were identified. The northern one was assumed to be the dominant branch of the SAF and the southern one was assumed to correspond to the PF. We find that the path of the SAF and PF closely track the underlying bathymetry. These jets narrow downstream of strong topographic features such as the Kerguelen Plateau, the Campbell Plateau, and the Pacific Antarctic Ridge, and they widen over the Bellingshausen Basin suggesting that topography may play a role in controlling the width of jets in these regions. Additionally, downstream of topographic obstructions, the SNRs of both the SAF and PF decrease, suggesting that lee effects cause enhanced eddy activity in these regions and/or that the jets become less stable.
Finally, we find that between 1993 and 2012, the ACC as a whole and on regional spatial scales has not exhibited significant meridional trends. The ACC was also shown to be relatively insensitive to climate mode-induced variability except in the East Pacific where SAM is correlated with changes in the basinaveraged SAF and PF position (a northward shift for positive SAM) and anticorrelated with the PF south of Australia (southward shift for positive SAM).