## Abstract

The mean state, variability, and extreme variability of the stratospheric polar vortices, with an emphasis on the Northern Hemisphere (NH) vortex, are examined using two-dimensional moment analysis and extreme value theory (EVT). The use of moments as an analysis tool gives rise to information about the vortex area, centroid latitude, aspect ratio, and kurtosis. The application of EVT to these moment-derived quantities allows the extreme variability of the vortex to be assessed. The data used for this study are 40-yr ECMWF Re-Analysis (ERA-40) potential vorticity fields on interpolated isentropic surfaces that range from 450 to 1450 K.

Analyses show that the most extreme vortex variability occurs most commonly in late January and early February, consistent with when most planetary wave driving from the troposphere is observed. Composites around sudden stratospheric warming (SSW) events reveal that the moment diagnostics evolve in statistically different ways between vortex splitting events and vortex displacement events, in contrast to the traditional diagnostics. Histograms of the vortex diagnostics on the 850-K (~10 hPa) surface over the 1958–2001 period are fitted with parametric distributions and show that SSW events constitute the majority of data in the tails of the distributions. The distribution of each diagnostic is computed on various surfaces throughout the depth of the stratosphere; it shows that in general the vortex becomes more circular with higher filamentation at the upper levels. The Northern and Southern Hemisphere (SH) vortices are also compared through the analysis of their respective vortex diagnostics, confirming that the SH vortex is less variable and lacks extreme events compared to the NH vortex. Finally, extreme value theory is used to statistically model the vortex diagnostics and make inferences about the underlying dynamics of the polar vortices.

## 1. Introduction

The winter stratosphere is dominated by a strong cyclonic vortex located in both the Northern Hemisphere (NH) and Southern Hemisphere (SH). Because of the abundance of landmass in the NH compared with the SH, more planetary waves are generated and the NH vortex is more variable than that of the SH (e.g., McIntyre and Palmer 1983). The edge of the vortices is defined by strong potential vorticity (PV) gradients that inhibit isentropic transport through the mechanism of “Rossby wave elasticity” (Juckes and McIntyre 1987). During extreme cases known as major sudden stratospheric warmings (SSWs) a distortion of the vortex can often lead to straining out and stirring of the vortex air mass into the background flow, resulting in deceleration of the cyclonic winds (e.g., Waugh et al. 1994). Such events are often defined as “major” when “the zonal mean zonal wind (ZMZW) at 10 hPa and 60°N reverse from westerly (eastward) to easterly (westward)” (Andrews et al. 1985) and “minor” if a significant weakening in the ZMZW is observed but with no reversal. During these extreme events the NH polar vortex has been shown to have a direct influence on tropospheric variability (e.g., Baldwin and Dunkerton 2001; Thompson et al. 2002; Bell et al. 2009). Furthermore, studies have shown that the winter stratosphere is in a highly disturbed state for around 10% of the time (e.g., Coughlin and Gray 2009). However, due to the lack of planetary waves in the SH, SSWs are almost nonexistent (e.g., Charlton et al. 2004; Esler et al. 2006) and most of the variability is observed in the upper stratosphere.

SSWs can be further subdivided into two types, with their own distinct structure and evolution of the vortex:

A

*vortex displacement*event is associated with anomalously high wavenumber-1 planetary wave activity entering the stratosphere and is characterized by a vortex with a comma-like shape that is shifting equatorward. Often this shifting occurs “top down” and the vortex has a baroclinic structure. Subsequently the Aleutian high, a weak anticyclone, encroaches over the pole and is especially dominant at lower levels.A

*vortex splitting*event is associated with anomalously high wavenumber-2 planetary wave activity entering the stratosphere. During such an event the vortex barotropically splits into two “daughter” vortices that tend to align along the 90°E–90°W axis, with one centered over Siberia and the other centered over northeastern Canada (Matthewman et al. 2009, hereafter M09).

Figure 1 shows these characteristics using the 850-K potential vorticity fields of the vortex during its (a) stable, (b) displaced, and (c) split states. Charlton and Polvani (2007) undertook a comprehensive classification of these types of SSW using the 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) and National Centers for Environmental Prediction (NCEP)–National Center for Atmospheric Research (NCAR) reanalysis datasets to determine the frequency of such events. They found that SSWs in the NH occurred on average twice every 3 years, with vortex splitting events occurring mainly in January and February and vortex displacements events occurring randomly throughout the winter period.

Because of the frequency of SSWs and their influence on the troposphere, many studies have characterized the structure and evolution of the vortex throughout the winter period (e.g., Manney et al. 1991; Limpasuvan et al. 2004). Most relevant for this study are the analyses using 2D-moment diagnostics undertaken by Waugh and Randel (1999, hereafter W99) and M09, which are themselves an adaptation of the methods used previously by others (e.g., Melander et al. 1986; Butchart and Remsberg 1986; Dritschel 1993).

W99 used 2D moments of the PV fields to examine the vertical and temporal evolution of the polar vortex during winter. M09 extended their moment analyses to include kurtosis and were able to determine whether or not the vortex was splitting during a SSW event.

The analysis undertaken in this study builds on previous work by further developing and refining the methods for calculating the moments, applying them to a longer period reanalysis dataset, and extending their interpretation to further examine behavior during extreme events. The characteristics of the NH moment diagnostics, especially during SSW events, are compared and contrasted with the more traditional diagnostics used to define SSWs, namely the zonal mean zonal wind at 60°N and the 60°–90°N polar cap temperatures. This comparison is carried out by compositing each of the vortex diagnostics around the central SSW dates defined by Charlton and Polvani (2007). The study is then extended by asking whether the distributions of the diagnostics can be modeled statistically, so that a set of simple distributions for the diagnostics can be identified from the reanalysis dataset that can then be used as a benchmark for the representation of the polar vortex in stratosphere-resolving models. This is done by fitting parametric distributions to each of the vortex diagnostics and testing their goodness of fit at various heights. A comparison of NH and SH hemisphere diagnostics is also included. The analysis suggests that standard parametric distribution fitting is adequate for the majority of the data, namely the days when the vortex is comparatively undisturbed. However, when the vortex is substantially disturbed, during SSW events, this approach works less well. The analysis is therefore extended to employ extreme value theory (EVT) in order to model and examine the behavior of the diagnostics at the tails of their distributions.

The layout of the paper is as follows. In section 2 the dataset and data manipulation is described. The methodology employed to derive the moments is described in section 3. The moments are then used to describe the climatological structure of the vortex during the whole winter period and during SSW events in section 4. Section 5 describes the parametric modeling of the distributions of the vortex diagnostics and section 6 describes the application of EVT to the tails of the distributions. Finally, in section 7, the results are summarized and discussed.

## 2. Data

This study uses ERA-40 reanalysis data obtained from the British Atmospheric Data Centre (BADC). Other similar studies have used NCEP data (e.g., Waugh 1997; W99); however, Charlton and Polvani (2007) undertook an analysis that compared SSWs from ERA-40 and NCEP–NCAR and concluded that there was very little difference in their representation. This held true throughout the 1958–2002 period that was studied. Additionally, Uppala et al. (2005) showed that the representation of SSWs in the NH winter was sufficiently accurate over the entire ERA-40 dataset (1958–2002). The SH, although less accurate, was also found to adequately represent the large-scale structure of the stratosphere. In particular, Uppala et al. (2005) commented on the similarities in the data between the SH 2002 major warming and the NH major warmings.

The ERA-40 dataset contains PV on pressure surfaces that extend to 1 hPa in the stratosphere. For the purpose of calculating moments PV is transformed onto potential temperature surfaces *θ* using the following:

where *T* and *p* are the atmospheric temperature and pressure respectively, *p*_{0} is the surface pressure of 1000 hPa, *R* is the molar gas constant, *c _{p}* is the specific heat capacity at constant pressure, and

*R*/

*C*is equal to for air (e.g., Ambaum 2010).

_{p}The PV is further interpolated onto standard theta surfaces that range from 450 to 1450 K at 200-K intervals using quadratic interpolation.

The change in coordinates from pressure to theta is necessary as PV is conserved under adiabatic flow on isentropic surfaces and can therefore be tracked on each surface. A further coordinate change from a polar to a Cartesian grid is also undertaken to simplify the calculation of moments (Waugh 1997); the results from these calculations are then mapped back to polar coordinates to give meaningful results. The transform to Cartesian coordinates is as follows:

where *λ* is the longitude and *φ* is the latitude in polar coordinates (Waugh 1997).

To further understand the structure and evolution of the vortex this study also analyses two traditional vortex diagnostics: the ZMZW at 60°N and the 60°–90°N cap temperature. These are determined on the following pressure levels, which roughly correspond to the theta surfaces used for the moments: 50, 20, 10, 5, 2, and 1 hPa.

## 3. Methodology

The structure and evolution of the polar vortices is analyzed using a comparison between traditional diagnostics that are commonly used in the literature and moment diagnostics (detailed in the previous section). The traditional diagnostics of ZMZW and 60°–90°N cap temperature can be trivially computed directly from the wind and temperature fields. However, the moment diagnostics require a more involved approach in their calculation, namely

the vortex edge must be defined and used to isolate PV in the vortex region from the rest of the NH;

the new PV field, which only represents the vortex, must be converted from polar to Cartesian coordinates for ease of calculating moments;

the moment algorithms must be run over the transformed PV field; and

the coordinate change is reversed for analysis of the results.

### a. Defining the vortex region

To calculate moments of regions that define the vortex, the vortex air must first be isolated from the nonvortex air. This ensures that only the structure of the vortex is considered in the moment analysis. To do this PV on the vortex boundary, *q _{b}* is calculated and subtracted from the entire PV field. All values that are negative (i.e., outside the vortex) are set to zero. This ensures that nonvortex air is not included in the moment calculations, and gives rise to the following modified PV field:

To define the modified PV field, *q _{b}* must first be calculated. This has been performed in different ways by W99 and M09. The vortex edge used in W99 was determined by calculating the mean PV at the location of the maximum meridional PV gradient over all winters in the NCEP dataset. This value was then used as the vortex edge for all days on each isentropic surface. M09 took a different approach to defining

*q*in order to show that their results were independent of the chosen vortex edge. They took

_{b}*q*as the average PV poleward of 45°N 9 days before the onset of a SSW, and then proceeded to use this value throughout the SSW life cycle. Although W99 and M09 both used fixed values for the vortex edge, this study chooses to calculate the vortex edge on a daily basis to better characterize the vortex region throughout the winter period. Applying this change to the methods of both M09 and W99 seem to give rise to unphysical results during extreme vortex variability. This is most likely due to diabatic processes changing the vortex edge during the breakup period of a SSW. Hence a different approach to calculating the vortex edge is required in this study.

_{b}In an unrelated study, Nash et al. (1996) developed an algorithm specifically to calculate the NH polar vortex boundary. This method is outlined below.

The PV field on an isentropic surface is transformed so that the highest PV isolines are centered around the pole and monotonically decrease toward the equator. The new PV field now has coordinates of longitude and “equivalent latitude” (McIntyre and Palmer 1984; Butchart and Remsberg 1986).

PV is plotted against equivalent latitude, resulting in a curve that is monotonic with steep regions at the subtropical surf-zone edge and the vortex edge, where the largest gradient represents the edge of the vortex, as detailed in McIntyre and Palmer (1983).

Often a case occurs in which two sharp gradients are observed. During such a case the vortex edge is further constrained by averaging the zonal wind along each PV isoline that is obtained from the plot described in (ii). In situations where two peaks in the PV gradient are observed there will often be only one peak in the average wind speed.

Nash et al. (1996) demonstrated that the true vortex edge could be determined from the peak in the PV gradient that corresponds most closely to the peak in the averaged wind, thereby obtaining

*q*._{b}

This method of calculating the vortex edge significantly increases computation time because of its high complexity; however, analysis of the PV fields suggests that the vortex edge using this technique is better defined during the breakdown of a SSW than previous methods. One additional adaptation that has been applied in this study over Nash et al. (1996) is that the PV fields at higher potential temperatures are smoothed using a 20° filter in the latitude and longitude directions before the edge is calculated. This smoothing acts to remove high levels of filamentation that occur in the upper stratosphere, thereby leaving mainly high PV that is associated with the main body of the vortex.^{1} It should be noted that in Nash et al. (1996) this problem was not addressed as they only calculated the vortex edge on the 450-K isentropic surface.

### b. Defining the moment equations

W99 and M09 used the vortex moment diagnostics to define an “equivalent ellipse” that was representative of the vortex on any given isentropic surface. To uniquely define an ellipse the following parameters are required: the ellipse area, centroid, and aspect ratio. Waugh (1997) showed that these quantities could be directly inferred from the moment equations. Equation (3) defines moments *μ _{ab}* of the modified PV field

*q**(

*x*,

*y*):

where *a* is the moment order in the *x* direction and *b* is the moment order in the *y* direction.

Equation (3) is adequate for deriving the zeroth- and first-order moments; however, a coordinate change must be applied for it to become meaningful for the higher-order moments. This change gives rise to the centralized moment equations, which calculate moments relative to the centroid of the modified PV field. This becomes useful for the higher-order moments in that, for example, the second-order moment will give the variance around the vortex center, rather than the variance around the original coordinate center. The centralized moment *J _{ab}* is defined as follows:

where and define the center of the vortex in the *x* and *y* directions, respectively.

The definitions (3) and (4) allow the vortex parameters to be calculated on a Cartesian grid. In this study the equations used to define the parameters are identical to those used in M09. These define the vortex area, centroid, and aspect ratio:

where *A* is the vortex area, is the vortex centroid location, and *r* is the vortex aspect ratio. Note that is converted back to polar coordinates to give the vortex centroid latitude, which shows how displaced from the pole the vortex is.

It should be noted that with the exception of the area, these are also the definitions used in W99. The area has been adapted by M09 to give an “objective area,” which gives both a measure of the vortex area and the vortex strength and can be regarded as a proxy for the vortex circulation. M09 used this measure as a substitute for the area diagnostic and demonstrated that it gave a good representation of the dynamical significance of the vortex on any given isentropic surface. We have also compared results from the objective area with standard area and found that they evolve in a similar fashion for a majority of the winter period, with the main difference occurring during extreme vortex variability where the objective area responds with a greater amplitude than the standard area. Hence in this study the objective area is chosen over the standard area.

### c. The kurtosis as a measure of SSW type

In addition to the ellipse parameters, M09 derives another parameter using the moment definitions. This parameter, the vortex excess kurtosis, is defined as follows:

where *r* is the vortex aspect ratio.

M09 uses this as a measure of the onset of a splitting event, during which the kurtosis drops below a threshold of −0.1. They further demonstrate that the vortex most commonly has a dumbbell-shaped structure at values of negative kurtosis, and that this more often than not is a precursor to a vortex splitting event. The threshold of −0.1 allowed M09 to calculate moments of each daughter vortex separately during a splitting event. Here, however, we do not fit diagnostics separately to each of the vortices, but rather examine how the moment diagnostic evolves over the entire NH regardless of the vortex state.

Expanding on M09, in the present study the kurtosis is used not only as a measure of splitting events, but also as a measure of displacement events. During a displacement event large amounts of PV often strip away from the main vortex, causing high filamentation. This tends to lead to a stronger than normal central core of the vortex surrounded by weaker PV, a structure that produces high values of kurtosis. Likewise, during a vortex splitting event the centroid location of the PV field will be between the two daughter vortices. This leads to a weak central core surrounded by higher PV, which gives rise to low or negative values of kurtosis.

The nonlinear nature of kurtosis is such that changes in the PV field will produce larger changes in the kurtosis than the other diagnostics; this comes from the fourth-order term in Eq. (8). Consequently one might expect the kurtosis to be highly dependent on the spatial resolution of the dataset. Figure 2 shows the evolution of the kurtosis over the winter of 1967/68 on the 850-K surface for ERA-40 (black), which can be considered representative of all winters. The gray lines show degraded versions of the ERA-40 PV field and indicate a good agreement with the original ERA-40 resolution. The kurtosis diagnostic seems to fall down when the grid resolution is degraded to 10° × 10° resolution (dashed line); this shows large deviations in the December–January transition, and in mid-January from the full-resolution field (black line). However, the spatial resolution of models tends not to be coarser than 5° × 5°, where Fig. 2 still shows good agreement with the full resolution of ERA-40.

One important feature of the winter presented in Fig. 2 is that it contains a wave-2 major SSW event that commences on 8 January 1968 (Charlton and Polvani 2007). This is clearly seen in the kurtosis evolution in Fig. 2, where the kurtosis drops below −0.1 at the beginning of January to define the split. However, the more dominant feature, which was not addressed in M09, is the large kurtosis after the splitting event, which is caused by one of the daughter vortices breaking up and mixing into the midlatitude air. Analysis of all large kurtosis events, as well as idealized cases, suggests that the following events in the PV fields lead to high positive kurtosis:

stripping of PV from the central core of the vortex;

the breakdown of a parent or daughter vortex, which in itself is associated with high levels of filamentation; and

a tightening of the central vortex core.

In reality this increase in vortex strength (iii) is never large enough to produce positive kurtosis values as high as features (i) and (ii) and often the large kurtosis values are a combination of two or more of the above points. The nature of kurtosis, and indeed all the vortex diagnostics, is studied in more detail over the subsequent sections.

## 4. The climatological structure of the stratospheric vortices

We begin by considering the mean state of the NH vortex throughout the winter period on the 850-K (~10 hPa) surface. This is followed by an analysis of how the vortex structure evolves vertically throughout the stratosphere. Finally we consider the structure and evolution of the vortex specifically during SSW events.

Figure 3 shows the climatologies of (a),(b) the traditional diagnostics of ZMZW at 60°N and 60°–90°N cap temperature on the 10-hPa surface and (c)–(f) the moment diagnostics of objective area, centroid latitude, aspect ratio, and kurtosis on the 850-K isentropic surface.

The ZMZW, cap temperature, and objective area all show strong seasonal cycles. Their seasonal cycles are a direct consequence of radiative heating from the sun, causing a higher equator–pole temperature gradient and hence strengthening the polar jet, especially during the December–January period. These three diagnostics are all strongly coupled through thermal wind balance, so a trend in one of the diagnostics will inevitably lead to trends in the other two. The climatology of the vortex centroid latitude indicates a smaller seasonal cycle that is not substantial compared to the variation of its interquartile range. Furthermore, the location of the vortex centroid is a function of the Rossby wave activity, and thus the apparent seasonal cycle could be an indication of nothing more than higher wave activity during January months. This is confirmed in the climatologies of aspect ratio and kurtosis, which both show more extreme positive values in the late January to early February period, during which both the aspect ratio and kurtosis indicate that the vortex is being stretched and higher filamentation is occurring. This is consistent with studies by, for instance, McIntyre (1982), who shows that planetary wave activity entering the stratosphere in this period is at a maximum for the year.

To determine how the vortex varies at all levels in the stratosphere, a February monthly mean of the six vortex diagnostics has been taken. Figure 4 shows the median (solid line) of (a),(b) the traditional diagnostics of ZMZW at 60°N and 60°–90°N cap temperature on surfaces from 50 to 1 hPa and (c)–(f) the moment diagnostics of objective area, centroid latitude, aspect ratio, and kurtosis on surfaces from 450 to 1450 K; the 0.25 and 0.75 quartiles are plotted either side as dashed lines.

Due to increased temperatures and wind speeds at high levels in the stratosphere, Fig. 4 clearly shows an increase in the mean February ZMZW, cap temperature, and objective area with altitude. The centroid latitude, however, does not show any large deviations with height, except at the lower levels of 450 and 550 K, which deviate from the rest of the vortex by a few degrees. A similar result was also obtained by W99 using a January averaged vortex.

The aspect ratio, perhaps surprisingly, shows a gradual decrease with height. Given the increased variability in the upper stratosphere, as confirmed by the high kurtosis values at upper levels, it might be expected that the vortex would have a higher aspect ratio at higher levels, but this is not seen. One hypothesis for this could be the presence of constant stationary planetary waves at the lower levels, causing an elongation of the vortex. It could also be due to the Aleutian high interacting with the vortex more at these lower levels.

A similar analysis to that shown in Fig. 4 was also performed on the vortex centroid longitude (not shown), a diagnostic that has been largely ignored in this study. The analysis showed no particular westward tilt with height, as M09 suggested might be the case, with a mean centroid longitude for the vortex between 550 and 1450 K of ~50°E. As with the centroid latitude, the centroid longitude deviated somewhat from the norm on the 450-K surface, suggesting that the weak PV gradients at this level are more dominated by Rossby waves than at other levels.

### Composite during SSWs

We continue this study by compositing the vortex diagnostics 20 days either side of the SSW central dates defined in Charlton and Polvani (2007), which are shown in Fig. 5 for vortex displacement events (solid black line) and vortex splitting events (dashed black line). The individual years are also plotted for splitting events (thin gray lines) and displacement events (thin black lines). The gray shaded regions show that the displacement and splitting events are statistically different at the 95% level according to a two-sample Kolmogorov–Smirnov test.^{2}

A comparison between the displaced and split composites for all diagnostics indicates that the moment diagnostics perform better than the traditional diagnostics at determining which type of SSW is occurring and this is confirmed by the statistical significance tests. It should be noted that Charlton and Polvani (2007) also found no statistical significance for similar traditional diagnostics when composited 20 days either side of the SSW onset, although they did not examine moment diagnostics in this fashion.

Perhaps the reason that the moment diagnostics are showing a higher contrast between splits and displacements is that they are essential measures of the geometry of the vortex, which is very different during the two types of event, whereas the traditional diagnostics use averaged quantities, which introduce unwanted information regarding the hemisphere and greatly reduces the amount of information regarding the vortex. Hence these traditional diagnostics may not be as sensitive to the vortex evolution as the moment diagnostics are. It should be noted that although in Charlton and Polvani (2007) splits and displacements were mostly calculated correctly using an algorithm based on the traditional diagnostics, the very nature of their work ensured that each SSW was either a split or a displacement. In reality SSWs exhibit wave 1 and wave 2 characteristics and therefore such a binary classification is not always appropriate. Evidence for this is in part seen from the individual winters plotted in Fig. 5. Using the composites of the moment diagnostics from Fig. 5, the following properties can be inferred:

The vortex centroid latitude is largely equatorward during the displaced composite (Fig. 5d, thick solid line). If we compare with the climatology of the vortex latitude (Fig. 3) we see that this equatorward shift is well outside the interquartile range at any point in the winter.

The vortex aspect ratio for the split vortex composite (Fig. 5e, dashed line) is significantly larger than the displaced composite. This in part comes from how we have defined our vortex. During a split vortex W99 only consider the larger of the two vortices, whereas M09 consider both vortices separately. Here we consider both vortices in one calculation; therefore, the moment diagnostics represent the shape of two separate vortices, producing a high aspect ratio. Using the aspect ratio in this way allows it to be a more useful diagnostic for indicating vortex splitting events than otherwise.

The kurtosis of the displaced composite (Fig. 5f, thick solid line) tends to high positive values, indicating strong filamentation is occurring, whereas it becomes negative during a splitting event (dashed line) with some filamentation that is observed around 10 days after the event. This delayed filamentation is a direct consequence of one of the daughter vortices breaking down via PV shedding (not shown). This analysis indicates that the kurtosis, in combination with the other diagnostics, can uniquely define the SSW type.

## 5. Modeling the distributions of vortex diagnostics

We continue this study by modeling each of the vortex diagnostics using parametric distributions. This has the advantage that the diagnostics can be compared, for example, on different heights or between the two hemispheres. In addition, it provides a framework for how the vortex may be represented in observational data and can therefore be used as an analysis tool to assess how well stratosphere-resolving models reproduce the vortex and its variability.

Figure 6 shows histograms over all winters [December–March (DJFM)] between 1958 and 2002 for the six vortex diagnostics. Note that the seasonal cycle has been removed from the ZMZW, cap temperature, and objective area. Note also that the centroid latitude data have been cubed to produce a Gaussian-shaped histogram because the original histogram did not accurately fit any known distribution. When interpreting the centroid latitude distribution, results are mapped back to normal latitudes for convenience. In addition to the histogram, the upper *x* axis indicates the average value of each diagnostic 3 days on either side of the central warming data for displacement events (black) and splitting events (dark gray).

The ZMZW and centroid latitude are fitted with a Gaussian distribution, which has the following probability density function (PDF) and parameters:

where *σ* is the standard deviation and *μ* is the mean of the data.

The remaining diagnostics of cap temperature, objective area, aspect ratio, and kurtosis are fitted with a generalized extreme value (GEV) distribution, which has the following PDF and parameters:

with

where *μ* is the location parameter that determines how far shifted along the *x* axis the distribution is, *σ* is the scale parameter that determines how spread out the data is, and *ξ* is the shape parameter that gives a measure of the skewness of the data (Cole 2001). These parameters are estimated using maximum likelihood methods to calculate the most probable value (see, e.g., Wilks 2006) and are detailed in Table 1 along with the standard errors in brackets.^{3}

To test how well the distribution fits the data, goodness-of-fit tests are employed. Here we choose to use the one-sample Kolmogorov–Smirnov test, which indicates that all the diagnostics are accurately modeled at the 90% confidence level. Note that other nonparametric tests are equally as valid for this analysis (e.g., Anderson–Darling and Cramér–von Mises tests; Choulakian and Stephens 2001). A qualitative approach to goodness-of-fit tests is achieved using quantile–quantile (Q-Q) plots, which provide a further, more in-depth test for each of the vortex diagnostics. These are shown beneath the corresponding histograms in Fig. 6. Q-Q plots essentially compare the quantiles of the diagnostics against the theoretical quantiles of the chosen distribution; accurate fits lie along the 1:1 line on Q-Q plots. The number of quantiles used to calculate the Q-Q plots is ~100, which is adequate to give a good overall representation of the datasets.

The Q-Q plots of ZMZW, cap temperature, objective area, and centroid latitude show that the data and distributions are well fitted. Small deviations are seen in the tails, but in all cases these deviations are less than 5% of the data. The ZMZW and cap temperature both show poor fits during the extreme phases that relate to SSWs: that is, <−40 m s^{−1} and >25 K, respectively. To accurately fit distributions to these extreme tails, a branch of statistics known as extreme value theory, which deals specifically with tail data, could be employed. This technique is used in the final analysis of the vortex diagnostics (section 6). A comparison between the histogram and Q-Q plot of the objective area demonstrates that data in the tails of the distributions in this case are also poorly fit. However, further analysis suggests that the data above 4 on the Q-Q plot all come from one event, a highly disturbed vortex over the winter of 1987/88 during which two SSWs occurred. After the recovery of such extreme events the vortex can often overintensify. This is due to the blocking of planetary waves at lower levels by the easterlies, allowing the vortex to tend toward radiative equilibrium. It is this process that gives rise to the high values of objective area seen on the Q-Q plot.

Moving on to the centroid latitude, the Q-Q plot seen in Fig. 6j, shows that the Gaussian distribution fits the data well except at both the upper and lower extremes. These extreme points only comprise around 5% of the total data, so the majority of the centroid latitude data are well fitted. The vortex displacement events (black dashes) seen on the histogram’s upper *x* axis also confirm that this diagnostic is especially useful in distinguishing these events from vortex splits. Here the displacement events are bunched together on the equatorward (lower) side of the distribution, whereas the splitting events are randomly distributed throughout.

Likewise, for the distribution of the aspect ratio (Fig. 6h) the splitting events are seen in the upper tail, confirming that it is a useful diagnostic for determining vortex splitting events. The kurtosis goes further in that there is a clear tendency for vortex displacement events to be in the upper half and vortex splitting events in the lower half, indicating its use as a diagnostic for identifying either type of SSW.

The Q-Q plots for both the aspect ratio and kurtosis (Figs. 6k,l) show good initial model-to-data fits, ranging from the lower tails to the beginning of the upper tails. The aspect ratio fits well to around 3.5, which, when compared with its histogram, encompasses 95% of its data. Likewise, for the kurtosis the distribution fits the data well up to a value of around 1, which encompasses >90% of the data. It should be noted that given the extreme nature of the positive kurtosis values, combined with a high number of data points, the Q-Q plot for the kurtosis indicates a bad model-to-data fit. However, in reality it is only the positive tail that has a bad fit. Furthermore, beyond a kurtosis of around 1 the Q-Q plot is very linear, indicating that the data fits a model similar to the GEV distribution multiplied by some scaling factor.

### a. Vertical comparison of the vortex

Having justified that the distributions detailed in Table 1 are indeed adequate to fit the bulk of each of the vortex diagnostics, the analysis is extended to fit distributions to the data on each of the isentropic surfaces from 450 to 1450 K at 100-K intervals for the moment diagnostics and from 50 to 1 hPa for the traditional diagnostics.

Each diagnostic is constrained to use the distribution named in Table 1 throughout the depth of the stratosphere. However, the parameters of each distribution are calculated directly from the data on each different surface and hence are allowed to vary; this is seen in Fig. 7.

The ZMZW and cap temperature both behave as expected in that the higher levels show more variability than the lower ones. As the seasonal cycles have been removed, absolute values cannot be directly obtained. However, as implied in section 4, the variability of the traditional diagnostics also increases with height; hence the upper stratosphere is more dynamically active.

This behavior is also observed in the objective area, where more variability is seen the farther up the diagnostic is in the stratosphere. It is also known that during SSW events the vortex often breaks down in a *top downward* fashion (e.g., Labitzke 1977). During these events planetary waves propagate along the vortex edge and break in the upper stratosphere, causing more variability in this region. If the planetary wave driving is strong enough this breaking will occur at subsequently lower altitudes, making variability in the strength of the vortex less common in the lower stratosphere. This behavior is clearly demonstrated in the distributions of objective area.

The centroid latitude, however, does not show any clear differences in the vortex structure throughout the stratosphere except on the 450-K surface, which suggests that the vortex is displaced the same amount from the pole throughout the stratosphere. At 450 K the vortex is weaker, and hence wave driving can often cause large displacements here. This result is confirmed in W99, who use a January average of the vortex centroid latitude and demonstrate that there is little deviation from the NH pole with height.

Perhaps surprisingly, the mean of the vortex aspect ratio (Fig. 7e) decreases with altitude, as does the variability. This is contrary to what might be expected, in that with more planetary wave activity in the upper stratosphere one might expect more warping of the vortex. However, this result is confirmed by studying PV and geopotential height fields of the vortex throughout the stratosphere. One explanation for this decrease in aspect ratio with height is the interaction of the Aleutian high, which has more influence at lower levels where it is encroaching onto the pole more (e.g., Harvey et al. 2002). In addition, radiative time scales in the lower stratosphere are ~30 days compared with ~5 days in the upper stratosphere (Newman and Rosenfield 1997); hence a vortex with a perturbation in its aspect ratio will take longer radiatively to return to its mean state.

Finally, the distributions of kurtosis show that filamentation of the vortex is occurring on average more at higher altitudes in the stratosphere, with the exception of the vortex on the 450-K surface, which shows a larger kurtosis. The higher filamentation that is occurring in the upper stratosphere results directly from more wave breaking at these levels, whereas the anomalous kurtosis distribution seen on the 450-K surface is most likely a consequence of the bottom of the vortex not always being well established, and therefore the breakdown and reforming of the vortex is more apparent.

### b. Comparisons with the Southern Hemisphere

This analysis is now further expanded to compare the NH and SH vortices on the 850-K and 10-hPa surfaces. The use of 2-sample Q-Q plots is now employed as an extra aid to this analysis. Here the Q-Q plot is used to compare two datasets instead of, in the case of the one-sample Q-Q plots, a dataset with a theoretical distribution. Figure 8 shows the distributions of the NH (solid) and SH (dashed) vortices for (a),(b) the traditional diagnostics of ZMZW at 60°N anomaly and 60°–90°N cap temperature anomaly on the 10-hPa surface (Figs. 8a,b) and the moment diagnostics of objective area anomaly, centroid latitude, aspect ratio, and kurtosis on the 850-K isentropic surface (Figs. 8c,g,h,i). Beneath each distribution the two-sample Q-Q plots comparing the Northern and Southern Hemisphere vortices are displayed. For the SH data the winter is defined as June–September (JJAS) and spans 1958–2001. Note that the final warming in the SH occurs around November–December (e.g., Haigh and Roscoe 2009) and therefore the JJAS months do not include periods where the vortex has permanently broken down.

The traditional diagnostics of ZMZW and cap temperature demonstrate more variability in the NH than the SH. The Q-Q plots of both diagnostics also indicate that it is the extreme negative winds and extreme positive temperatures that are causing the largest differences between the two hemispheres. This is due to the occurrence of SSW events in the NH, which are largely nonexistent in the SH.

The distributions of the vortex objective area also show a bias toward the NH being weaker. This is confirmed in the aspect ratio Q-Q plot and caused by the presence of more planetary waves in the NH stratosphere. It should be noted that because the data have had the seasonal cycle removed, its mean is by definition 0; hence the (*x* = 0, *y* = 0) point on the Q-Q plot must lie on the 1:1 line.

From the distributions of the centroid latitude it is immediately clear that the polar vortex in the SH is more poleward than that of the NH, with a mean centroid latitude of 85°S in the SH compared to the NH mean latitude of 77°N. This is confirmed in both the centroid latitude Q-Q plot and in the study undertaken by W99, who found that the vortex in the SH was located significantly closer to the pole than in the NH. The mean centroid latitude of the NH vortex in this study and that of W99 was 77° and 76°N, respectively, showing good agreement. Similarly, the mean centroid latitude of the SH vortex from this study was 85°S, compared with 86°S in W99. Perhaps the reason for the slight disagreement is the addition of September in our choice of the SH winter, as well as the number of years in each analysis. The analysis used in W99 extended from 1979 to 1998 using the June–August (JJA) months, where as in this study the JJAS months from 1958–2001 are used. Note that the only known SH warming, which occurred in September 2002, is not included because ERA-40 does not extend this far.

The distributions of both the aspect ratio and kurtosis show consistently higher values and more variability in the data for the NH than the SH. Furthermore the Q-Q plots of both diagnostics indicate that the more extreme the data are in the upper tails, the bigger the differences between the two hemispheres. This is seen in the rate of the linear increase in the Q-Q points with respect to the 1:1 line in both the aspect ratio and kurtosis. In comparison with W99, who obtained a mean aspect ratio in the SH of 1.2, the mean aspect ratio obtained here is 1.3, again indicating that the period used in this study of JJAS is more variable than JJA.

## 6. Extreme value modeling

Much of the analysis presented in this study has indicated that the modeling of the extreme data is badly represented by the distributions (see in particular Fig. 6). In addition it was shown that the extreme values often correspond to SSW events and therefore are of greater importance to tropospheric and stratospheric teleconnections. We therefore proceed to model the extreme values of the diagnostics in order to better understand these events. This type of analysis is undertaken using a branch of statistics known as extreme value theory, which deals specifically with fitting distributions to data over a certain threshold.^{4} One significant departure from the analysis of the complete datasets is that data here are declustered, allowing each data point to be independent. Each cluster of a time series is defined as consecutive exceedences over the chosen threshold *x*_{th}. The maximum data value within each cluster is taken as a single independent data point to represent the cluster as a whole. This technique prevents the modeling of extremes from being dominated by single events and instead allows the distribution to be representative of all events. The method of choosing such a distribution to fit the declustered data is similar to the previous section (see also Cole 2001). For this study the generalized Pareto (GP) distribution is used for all vortex diagnostics. This distribution has the following PDF:

where the distribution of data *x* is defined by two parameters: *ξ*, the shape parameter, and *σ*, the scale parameter. (For more discussion of these parameters see the previous section.) Specifically with regard to the GP distribution the shape parameter *ξ* can also be used to determine whether or not the dataset has a physical limit. If *ξ* < 0 then the GP distribution is unbounded for the upper tail; therefore data can tend to infinity. If *ξ* ≤ 0, then the GP distribution has an upper limit. The maximum value of the upper tail is then predicted using the following relation (Cole 2001):

Figure 9 (first column) shows the analysis of fitting the GP distribution to the vortex diagnostics, (second column) the corresponding Q-Q plots, and (third column) the return level plots, where the return level plot is the average frequency for a given value to recur. The parameters of the GP distribution, threshold, and 10-yr return level are all summarized in Table 2. Note that the 60°–90°N cap temperature diagnostic has been dropped from this analysis as it is highly correlated to the ZMZW at 60°N, and the previous analysis has shown that it provides little or no extra information regarding the structure of the vortex.

The GP distribution fits the ZMZW anomaly at 60°N well. A minor misfit is apparent in both the density and Q-Q plots between −35 and 40 m s^{−1}; however, the broad structure of the ZMZW anomaly is captured well by the model. Both the negative shape parameter and curvature of the return level plot suggest that the ZMZW anomaly has a physical limit. This limit can be determined by Eq. (13) and is calculated to be −68.9 m s^{−1}. A comparison with the ZMZW at 60°N for the winters outside the ERA-40 time frame (2002/03–2009/10) using data from the UK Met Office shows that this limit has not been exceeded, even during the comparatively disturbed winters during this recent period.

Moving on to the moment diagnostics, the GP distribution fit to the vortex objective area (Fig. 9, second row) is not as good as the remaining diagnostics. This mainly comes from the data at ~−2 × 10^{18} PVU m^{2}; the remaining data are reasonably well modeled by the distribution. The shape parameter for this diagnostic is not negative, consistent with the fact that the vortex cannot have an area less than zero and hence has a physical limit. However, the shape parameter of 0.08 is very small compared to the standard error of 0.19. Therefore the parameter could be positive or negative and gives more credibility to the GP distribution fit in this case.

The centroid latitude is again well modeled using the GP distribution, although the Q-Q plot suggests that values below 60°N do not fit the model so well. For this analysis 18% of the centroid latitude data are above the threshold, which is higher than might be expected for extreme events. However, the average cluster size for the centroid latitude is 10 events; therefore, the number of exceedences used to model the GP distribution is similar to those of the other diagnostics. A large grouping of the data is observed between 72° and 70°N. Analysis of the original clusters for this grouping indicate that the centroid of the vortex is often displaced to these latitudes for shorter periods (2–3 days) than the remaining extreme data. This suggests that the vortex is often displaced but not destroyed, consistent with theories on Rossby elasticity (e.g., Juckes and McIntyre 1987). Analysis of the return level plot suggests that once in every 10-yr period the centroid latitude is predicted to reach ~57°N, with a minimum predicted possible latitude within the climate system of ~13°N [using Eq. (13)]. This suggests a physical limit to how far the planetary waves can disturb and displace the vortex and although 13°N seems unrealistic, the prediction is most likely referring to a highly weakened vortex that has undergone major deformation. Additionally, the return period for a latitude of 13°N is well over a thousand years and therefore is an extremely rare event.

For the aspect ratio all three diagnostic plots in Fig. 9 show an exceptional fit between the GP distribution and the data. In particular the values above 4 on the Q-Q plot normally deviate greatly from the 1:1 line for even small discrepancies in the data. The curvature of the return level prediction, coupled with the positive shape parameter estimate for the aspect ratio, indicates that very large aspect ratio events can be reached and are not constrained beyond a certain value. Finally, the density plot indicates that the GP distribution fits the data well throughout the entire tail. Therefore in the case of the aspect ratio the GP distribution proves to give an excellent overall representation of the data, and hence provides an accurate benchmark as to how this diagnostic should be represented in other datasets.

Likewise, the three diagnostic plots for the vortex kurtosis (Fig. 9, bottom row) show high levels of agreement between the model and data, with the exception of the most extreme event, caused by the high levels of nonlinearity of the diagnostic. The modeling of the tail using the GP distribution shows great improvement over that of the GEV analysis (see previous section), which could not model the data above a kurtosis of 2. The use of the GP distribution, however, models the tail well up to a kurtosis of ~10. In addition, the kurtosis only exceeds its threshold on average for 3 days at a time, implying there are more short-lived kurtosis events compared to the other vortex diagnostics. This ties in well with the nonlinear nature of the diagnostic and suggests that PV is stripped from the main body of the vortex in short-lived, frequent events.

Approaching the extreme data of the vortex diagnostics from a SSW framework, the maximum aspect ratio observed for the vortex is 5.9 and occurs very closely to a vortex split on 22 February 1979. Additionally the centroid latitude of the vortex reached 51°N ~12 days after the onset of a vortex split, which occurred on 25 February 1999. During this event one of the daughter vortices is mixed into the midlatitude PV and leaves one highly displaced vortex remaining. It should be noted that for the majority of splitting cases the stronger vortex will often propagate poleward before the other vortex is mixed out. Consequently the vortex split over February 1979 was truly a rare event. It is interesting to note that many of the points classed as extreme data in the vortex diagnostics are not directly related to major SSW events. One example occurs on 2 February 1995, which shows the centroid latitude of the vortex located at 54°N (on the 850-K surface). This is comparable to the largest known displacement mentioned earlier and is of particular interest because the scientific community often considers the 1990s to be a largely undisturbed period for the polar stratosphere (e.g., Charlton and Polvani 2007). Observations of the remaining vortex diagnostics around this period suggest that associated filamentation occurs, as well as an elongation of the vortex ~10 days prior to the event. This particular event is however by no means unique. Many other large events not corresponding to major SSWs have been identified in the extreme data. Such analysis calls into question the use of the WMO definition to identify SSW events and should prompt further investigation into the characterization of the polar vortices.

For all the diagnostics the extreme data comprise on average ~100 events, corresponding to approximately two events per year over the ERA-40 dataset. This frequency is 3 times higher than that of SSW events predicted by Charlton and Polvani (2007). Naturally the question that arises is whether or not these extremes in the moment diagnostics can be correlated to changes in tropospheric climate, much in the same fashion as Baldwin and Dunkerton (2001) showed that strong/weak vortex events had an influence on surface climate. Although this question is beyond the scope of this current study, we hope to address it in future analyses.

## 7. Summary and discussion

We have compared six vortex diagnostics of ZMZW at 60°N, 60°–90°N cap temperature, objective area, centroid latitude, aspect ratio, and kurtosis to provide a comprehensive analysis of the structure and evolution of the vortex throughout the winter (DJFM in the NH and JJAS in the SH). Our calculation of moments differs from those of W99 and M09 mainly in our definition of the vortex edge, where we use a method laid out in Nash et al. (1996) to calculate the PV on the vortex edge on a day-by-day basis. This definition has led to moment calculations being robust enough to determine the diagnostics during many different states of the vortex, including the most extreme cases of SSW events. Using this method the following conclusions have been obtained:

The moment diagnostics of objective area, centroid latitude, aspect ratio, and kurtosis behave in different ways for vortex displacement events compared to vortex splitting events and can therefore be used to distinguish between them, in contrast to the more traditional diagnostics based on zonal wind and polar temperatures. This is especially true around the central data of an SSW [defined in Charlton and Polvani (2007)] as confirmed at the 95% confidence level using a Kolmogorov–Smirnov test. This implies that the moment diagnostics are a better measure of vortex variability than the traditional methods, although we recognize that the combination of both sets of diagnostics is preferable.

The centroid latitude can be used as a diagnostic that determines vortex displacement events, whereas the aspect ratio determines vortex splitting events. In addition the kurtosis can be used to pick out either type of event, where a large positive kurtosis indicates filamentation and a negative kurtosis indicates a splitting of the vortex.

The variability of the vortex as diagnosed by its zonal wind speeds, temperature, area, and strength increases with height as does the frequency of filamentation. The centroid latitude and longitude of the vortex, however, show no real deviations from the pole with height, with the exception of the less well-defined vortex on the 450-K level. Finally we note that the vortex becomes more circular at higher levels.

The vortex diagnostics confirm the reduced variability and lack of extreme events in the SH compared with the NH vortex. In the most extreme events in the NH the vortex can be centered around 50°N compared with 70°S in the SH. Likewise the aspect ratio of the vortex in the NH can reach between 5–6, compared with a value closer to 2.5–3 in the SH. These results also act to confirm those from W99, who found the mean latitudes over the winter period of 76°N and 86°S in the NH and SH, respectively, whereas in this study the corresponding values are 77°N and 85°S. W99 also calculated the mean of the aspect ratio over the winter periods and obtained 1.7 and 1.2 for the NH and SH, respectively, again close to the values of 1.7 and 1.3 that were obtained here.

The extreme value analysis of the vortex diagnostics in the midstratosphere provided a further, more accurate benchmark to modeling extreme vortex events and should be used in combination with the analysis using the Gaussian and GEV distributions. Such extreme events were found to often relate to major and minor SSWs. The fitting of a GP distribution to the extreme data of the vortex diagnostics also predicted that an absolute limit of the centroid latitude of the vortex was observed at 13°N, but that large displacements to ~57°N could be achieved once every 10 years. Associated with such events, ZMZW anomalies were predicted to have a maximum constraint of ~70 m s

^{−1}, consistent with analysis using data from the Met Office over periods outside the ERA-40 time frame.

One of the major conclusions drawn from this study is that traditional definitions of SSW events often miss extreme behavior of the vortex. The analysis of the distributions of each of the moment diagnostics has shown that so-called minor SSW events often have a larger influence on the geometry of the vortex than major SSW events. Likewise, the subdivision of splitting and displacement events is not always as well defined as the names suggest. Often a vortex displacing event will split into two daughter vortices before the vortex returns to its stable state, and vice versa. It is therefore not only the use of such traditional diagnostics that often leads to poor characterization of the vortex, but the very binary definition of SSW events that is used in the mainstream literature.

### Further analysis

The typical definition of a SSW, that the ZMZWs at 60°N and 10 hPa must reverse, can often result in disregarding large variability of the vortex. In many cases the ZMZWs at 60°N and 10 hPa will decrease to values just above the zero line and hence be classed as minor warmings. However, in reality they are often very similar to major warming events. Coughlin and Gray (2009) used the *k*-means clustering technique as a method of distinguishing between disturbed and undisturbed vortex states. However, using the traditional diagnostics of zonal mean winds and temperatures they found that the stratosphere was naturally split into only two states using this technique (i.e., they could not further subdivide the disturbed vortex state into splits and displacements). However, Hannachi et al. (2011) have used a similar clustering technique applied to the moment diagnostics instead of the traditional diagnostics and revealed three states corresponding to two types of SSWs and the stable vortex, thus showing the added benefits of the moment diagnostics.

Given the nature of these diagnostics, often large vortex variability is observed without a ZMZW reversal. It is the structure and evolution of the vortex that the moment diagnostics measure. During external forcings, for example El Niño, observational studies have shown that the vortex is in a disturbed state (e.g., Garcia-Herrera et al. 2006; Manzini et al. 2006). This is confirmed extensively in modeling studies (e.g., Sassi et al. 2004; Bell et al. 2009). However, the structure of the vortex during such events has not been examined closely. Using moment-based diagnostics not only during El Niño events but also during other external forcings such as the quasi-biennial oscillation, solar cycle, and large volcanic eruptions would provide a more detailed response to such forcings and hence a better understanding of vortex teleconnections.

One of the aims of this study was to provide a statistical representation of the vortex throughout the ERA-40 dataset. This can be seen in Fig. 6, details of which are noted in Table 1. The analysis here has shown that each diagnostic can be modeled parametrically, and that these distributions can be applied to data on a range of surfaces (e.g., hemispheric or isentropic). This versatility of the moment diagnostics using the methods laid out in this paper should mean they can also be applied to model output and hence test how well the vortex is represented. This could prove to be a useful analysis tool when attempting to validate stratosphere resolving models. It could also be expanded on by modeling the tails, or “extremes,” of the distributions separately to better characterize the extreme vortex variability in the models.

## Acknowledgments

DMM is funded by the Natural Environment Research Council (NERC). LJG is funded by the NERC National Centre for Atmospheric Science (NCAS). We thank N. J. Matthewman, J. G. Esler, A. Hannachi, C. Ferro, and C. Bell for earlier comments on these results. We would also like to thank the two anonymous reviewers for there helpful comments and suggestions regarding this paper.

## REFERENCES

## Footnotes

^{*}

Current affiliation: Atmospheric, Oceanic and Planetary Physics, Department of Physics, University of Oxford, Oxford, United Kingdom.

^{1}

The smoothed PV field is only used in the calculation of the vortex edge whereas moments are applied to the full PV field. In addition, smoothing filters were also applied at lower vertical levels to ascertain the sensitivity of smoothing to the vortex edge calculation. The results showed very little change between the smoothed and nonsmoothed PV fields. In principle a 20° filter could have been applied to PV on all levels, although this would have increased computational time.

^{2}

The Kolmogorov–Smirnov test compares how well two datasets agree with each other. Unlike other significance tests (e.g., the *t* test), it does not assume that the underlying data fit a known distribution.

^{3}

A standard error is essentially a standard deviation of the sampling distributions, which in turn is built from all possible values of the parameter estimate.

^{4}

The process of choosing a threshold is not straightforward. The reader is referred to Cole (2001) for more details on how this is performed.