CycloTRACK (v1.0) – tracking winter extratropical cyclones based on relative vorticity: sensitivity to data filtering and other relevant parameters

In this study we present a new cyclone identification and tracking algorithm, cycloTRACK. The algorithm describes an iterative process. At each time step it identifies all potential cyclone centers, defined as relative vorticity maxima embedded in smoothed enclosed contours of at least 3× 10−5 s−1 at the atmospheric level of 850 hPa. Next, the algorithm finds all the potential cyclone paths by linking the cyclone centers at consecutive time steps and selects the most probable track based on the minimization of a cost function. The cost function is based on the average differences of relative vorticity between consecutive track points, weighted by their distance. Last, for each cyclone, the algorithm identifies “an effective area” for which different physical diagnostics are measured, such as the minimum sea level pressure and the maximum wind speed. The algorithm was applied to the ERA-Interim reanalyses for tracking the Northern Hemisphere extratropical cyclones of winters from 1989 until 2009, and we assessed its sensitivity for the several free parameters used to perform the tracking.

The typical methods for cyclone detection and tracking follow a two-step approach.First, they identify the locations of cyclone centers at each time step and then the cyclones paths are extracted by connecting the centers of consecutive time steps.During the identification step, several constraints are applied in order to control the range and the number of the identified features.For example, in some studies, the definition of the location of a cyclone implies three constraints on the fields of mean sea level pressure: (1) the representative grid point of the data field has to have the minimum value among the adjacent grid points, (2) the minimum value has to be inferior to a threshold value, and (3) the field gradient has to be superior to a threshold value (e.g., Murray and Simmonds, 1991;Blender and Schubert, 1997;Nissen et al., 2010).However, there is a trade-off.The application of such "strict" constraints on pressure gradients may lead to tracking cyclones only during their mature stage and, furthermore, weak cyclones may be left undetected.
An efficient cyclone-tracking algorithm has to be able to decide whether the identified features have moved over time or whether they have ceased to exist.An extra challenge is that cyclones can split or even merge with other cyclones.Also, there might exist more than one candidate location to E. Flaounas et al.: CycloTRACK (v1.0) be considered as the next cyclone center.This is the case of noisy fields, where a significant number of grid points, typically located close to each, can be considered as cyclone centers.The algorithm has to determine automatically which of the candidate features should be tracked and which should be neglected.Many methods apply an across-time "nearestneighborhood" approach where tracks are built by connecting the identified cyclone centers of one time step with the nearest neighbors of the following time step (Blender et al., 1997;Serreze et al., 1997;Trigo et al., 1999).Other studies use more complex tracking algorithms based on displacement speed (e.g., Murray and Simmonds, 1991;Wernli et al., 2006;Davis et al., 2008;Campins et al., 2011;Hanley and Caballero, 2012).These algorithms make a "guess" as to the next step location of the cyclone and choose the nearest feature detected at that potential location.Finally, Inatsu (2009) presented an algorithm where cyclones are identified as areas of enclosed grid points that satisfy a certain condition on filtered meridional wind at 850 hPa.Then, cyclone tracking is performed by connecting the cyclone areas that overlap in consecutive time steps Hodges (1999) proposed post-processing of the identified features and their tracks.His method starts by building all potential tracks using the nearest-neighborhood approach.Then, the identified tracks exchange track points between them until an overall cost function is minimized.This cost function is a measure of the smoothness of the whole set of tracks.Hanley and Caballero (2012) also applied a postprocessing approach in order to cover the cases where cyclones that have more than one center undergo any merging or splitting process.Raible et al. (2008) were the first to compare the performance of three different tracking methods, applied to extratropical cyclones.Their results converged on the interannual variability of cyclone occurrences; however, they differed on the cyclone number trends and track densities.Recently the IMILAST (Intercomparison of mid latitude storm diagnostics) project presented a comparison of the performance of 15 different algorithms which have been applied for the tracking of extratropical cyclones during the cold season of 21 years over the entire planet (Neu et al., 2013).Neu et al. (2013) found that the number of tracks, the lifetime and the intensity of cyclones may vary significantly depending on the algorithm used.This apparent disagreement of the algorithms can be easily explained by the fact that there is no commonly accepted definition of a cyclone.Consequently, each algorithm applies different constraints and/or different fields.In this sense, one of the main results of Neu et al. (2013) is that no algorithm can be considered to be "superior" or more "correct" than the others, since they use different definitions of a cyclone.However, it is noticeable that even algorithms with similar configurations may present a divergence of results.Recently, Ulbrich et al. (2013) treated the evolution of extratropical cyclone tracks in the context of a changing climate using a multi-tracking-method approach.It was shown that, despite the differences on the tracks number, the algorithms present comparable results on the cyclones climatology future trends.This finding confirms that there is a common robust behavior independent of the different algorithms configurations and modeling constraints.
In this study, our principal motivation was to design an algorithm that is able to provide qualitative characteristics of the features like splitting, merging, wind speed, associated rainfall and minimum pressure in parallel with the tracking.A new aspect of the proposed approach is that cyclonic features are tracked based on their physical properties by assuring a gradual evolution of the cyclone relative vorticity, not on their displacement.The use of relative vorticity presents some advantages compared to the use of geopotential height or mean sea level pressure.It is a high-frequency variable, representative of local scales that -presumably -permit cyclone tracking at an early stage, i.e., during its initial perturbation and before it is characterized by closed pressure contours (Sinclair, 1994(Sinclair, , 1997;;Hodges, 1999;Inatsu 2009;Kew et al., 2010).This can be an advantage in cases where cyclones intensity increases significantly within 24 h, i.e., explosive cyclogenesis (e.g., Sanders and Gyakum, 1980;Trigo, 2006;Lagouvardos et al., 2007).The disadvantages of relative vorticity is that it is a wind-based field, sensitive to the horizontal resolution of the data set, and that local maxima might not correspond to wind vortices but to other features such as an abrupt wind turning.
To deal with the spatial noise of relative vorticity, we smooth the input fields.The smoothing operation partly counteracts the property of relative vorticity to detect cyclones in their early stage; however our algorithm has a high degree of flexibility that permits the tracking of perturbations that did not evolve into strong cyclones.A similar approach has been used before in other studies for capturing weak cyclonic features (e.g., Murray and Simmonds, 1991;Pinto et al., 2005), but in our case it provides an additional value for optimizing the algorithm and determining the cyclones that are not sensitive to filtering.The application and assessment of our method is performed in line with the efforts of the IMILAST project.We used the same time periods and input data sets in order to make our results comparable with those of the aforementioned project.
In Sect. 2 we describe cycloTRACK in detail, our cyclone detection and tracking algorithm.In Sect. 3 we present the results of several sensitivity tests of our method, applied to the ERA-Interim (ERA-I) data set for the winters (December-January-February) of the period 1989-2009. Finally, Sect. 4 presents the conclusions and our prospects for future research.

Identification and tracking algorithm
In this section we present our algorithm, cycloTRACK, and its application to the vorticity fields at 850 hPa level within  et al., 2008).The algorithm is composed of two independent steps.In the first step, we identify all cyclonic features for all time steps of a data set, and in the second step we build the cyclone tracks.

Step 1: identifying cyclones and quantifying their characteristics
The first step of the algorithm is devoted to the identification of the cyclones and the quantification of their characteristics.Initially, the algorithm identifies all cyclonic features, or more precisely all cyclonic circulations.Then, for each cyclonic circulation, the algorithm identifies the representative centers which will be treated as distinct cyclones.Finally, for each center, the algorithm quantifies its characteristics, such as the maximum relative vorticity, the maximum wind speed and the minimum sea level pressure.

Identification of cyclonic circulations
To identify cyclonic circulations, the vorticity field is smoothed by applying a spatial filter.In previous studies the vorticity field has been smoothed by a variety of filtering operations such as B-spline techniques (Hodges, 1995), time band-pass filtering (Hoskins and Hodges, 2002;Inatsu, 2009) and 1-2-1 filters (Satake et al., 2013).In this study, we use the simple method of a 1-1-1 spatial filter, which efficiently smoothes the orographic or coastal vorticity maxima and the gradients of relative vorticity fields.The latter also helps the algorithm to reject local vorticity maxima that are nested within noisy field gradients -especially when considering very high resolution data sets.The smoothing operation on the relative vorticity field is performed at each grid point separately by multiplying the sum of all its adjacent grid points within distance X by 1/(2X + 1) 2 .For instance, at any grid point (a, b) the smoothed relative vorticity (RV) is given by RV a,b = 1 As X increases, the smoothing operation on the relative vorticity field becomes stronger.Finally, we apply a threshold and we retain only the grid points exceeding this threshold.
Figure 1 shows the raw relative vorticity fields and the corresponding fields after the application of three different filters with 2X + 1 equal to 3, 5 and 7.The relative vorticity fields are derived from ERA-I and are centered over Europe at 00:00 UTC, 3 December 1999, featuring the storm "Anatol" over Denmark as the strongest detected cyclone.In all panels of Fig. 1, the threshold is set to 3 × 10 −5 s −1 .As the applied filter becomes stronger, the relative vorticity values become weaker.Small vorticity features tend to be suppressed; however the structure and location of the vorticity maxima of the strongest features, such as Anatol, are not altered among the different filter operations.We used filtering for smoothing the values within a cyclonic circulation.Thus, the filtering matrix should not be much larger than the length scale of a cyclone.In this sense, a 7 × 7 grid point filter for ERA-I means that relative vorticity is smoothed over a 10.5 • × 10.5 • region, which is a rather large area.As shown in Fig. 1a and b, each cyclonic circulation might correspond to a unique cyclone or to a large complex of cyclonic centers comprised by more than one local maximum.The 3 × 10 −5 s −1 threshold applied to the ERA-I data set (1.5 • × 1.5 • resolution) has been found adequate for describing cyclones, even at their initial stage, for all three filtering sensitivity tests.In this step, the algorithm identifies and labels, with a number, all cyclonic circulations which are defined as the areas composed of neighboring grid points of values exceeding the 3 × 10 −5 s −1 threshold.The threshold value allows for the algorithm to be tuned for detecting cyclones in coarse resolution data sets (e.g., 1.5 • × 1.5 • , as in ERA-I used here) or in high-resolution data sets (e.g., 20 km regional climate runs).Applying a threshold is a convenient approach for adjusting the filtering strength.Alternatively, we could keep the filtering strength constant and allow the threshold value to vary.However, it is only by varying the filtering strength that the vorticity field may be smoothed within the characteristic length scale of cyclones.Similar approaches for identifying a feature through an enclosed area have been used before for cyclones (e.g., Hodges, 1999;Wernli and Schwierz, 2006;Inatsu, 2009;Flaounas et al., 2013) and for MCSs (e.g., Machado et al., 1998).

Identification of cyclonic centers
Careful inspection of Fig. 1b, c and d reveals that the cyclonic circulations do not correspond to the same cyclone.For this reason, each labeled cyclonic circulation is further treated in order to locate all embedded local vorticity maxima.These local maxima will be labeled, and eventually they will be treated as centers of unique cyclones.The term "centers of unique cyclones" has no physical basis but it is conveniently used here in order to describe the grid points that present local maxima of relative vorticity and that we follow in time in order to construct cyclones tracks.In this sense we need to provide the algorithm with a representative cyclone center even though the cyclone structure might be very complex, with more than one vorticity maximum (especially in very high resolution data sets).To deal with this issue we (1) filter the data, smoothing the noisy gradients (already performed in the previous step); (2)define the local maximum as the maximum value of the central grid point among its eight surrounding grid points; and (3) consider that, between two centers, there is a relative vorticity difference greater than a threshold value (in this case set equal to 3 × 10 −5 s −1 ) which is applied to define the cyclonic circulations.The last criterion prohibits weak cyclonic circulations, i.e., cyclones of relative vorticity close to the threshold value, in order to present multiple centers.

Quantifying cyclone characteristics
Once all cyclones have been identified, we determine an "effective area" for each cyclone.This area is a circular disk centered at the cyclone vorticity maximum.The disk radius grows gradually until (1) all grid points included in the disk have a vorticity average inferior to a threshold value, (2) the radius reaches a pre-defined maximum length or (3) a relative vorticity value greater than that of the cyclonic center is found within the area.According to this empirical method, the strong or the large and weak cyclones tend to produce large effective areas.The third criterion favors stronger cyclones spreading their area independently of the presence of other weaker ones in their region, while it prevents the weaker cyclones from sharing the same area with stronger cyclones.In Flaounas et al. (2013), the cyclone area was defined by the cyclone-enclosed contour as defined by the applied threshold value (see their appendix figure).However, such an enclosed area might not capture grid points that present relative vorticity values lower than the applied threshold.In Lim and Simmonds (2007), the cyclone area was defined by a representative circular disk of radius equal to the average distance between the cyclone center and the enclosing zero contour of the mean sea level pressure Laplacian.The circular disk seemed to be the best choice for our algorithm in order to capture the areas affected by a cyclonic vortex.Irregular shapes may also be considered, such as instance enclosed contours of pressure (Wernli and Schwierz, 2006;Hanley and Caballero, 2012) or relative vorticity (Flaounas et al., 2013).
Once the effective area is defined, our algorithm computes the physical properties of the cyclone within it.Figure 2 shows, as an example, the effective area and the detected minimum sea level pressure and maximum 10 m wind of Anatol at the same time as in Fig. 1b.

Step 2: tracking cyclones
In this step we combine the cyclone centers into a track.First, the algorithm sorts the identified cyclones based on their relative vorticity value, from the strongest (i.e., the one with the highest relative vorticity value) to the weakest.Then, it starts from the first cyclone and searches forward and backward in time for all its possible tracks.More precisely, the algorithm constructs all possible cyclone tracks that share the same maximum vorticity track point.Once all possible tracks are constructed, the algorithm chooses the track that presents the most "natural evolution" of relative vorticity, i.e., the track which presents the smallest differences of relative vorticity in consecutive points, weighted by the distance between the track point locations.Figure 3a illustrates an idealized experiment, presenting the locations of all identified cyclones in a four-time-step data set.Six cyclones are identified: one in the first time step, one in the second time step, and two each for time steps three and four.The tracking process begins from the strongest cyclone (i.e., the cyclone 2(12)) and constructs all possible tracks by iterating forward and backward in time with all other features.Figure 3b shows that the first cyclone may undertake four possible tracks; however it is obvious that tracks 1(9), 2(12), 3(10) and 4(8) present the most "natural evolution", since maximum relative vorticity presents the smallest difference from one time step to the next.The algorithm saves this track and deletes the locations of the cyclones used from the data set.Then, a new iteration begins in which the algorithm starts from the cyclone with the highest vorticity and a new track is eventually constructed (Fig. 3c).We found that starting the tracks from the cyclone's mature state is more efficient for the construction of the first steps of the tracks.Indeed, for most cases in the previous and in the following time step of the cyclone with the highest vorticity state, there is only one strong cyclone to act as a candidate for continuing the tracks.
The practice of cost function minimization has been used in relevant literature of tracking algorithms.Hodges (1995) built the feature tracks by minimizing the cost function of the feature's track smoothness, while Hewson and Titley (2010) built the feature tracks by applying a likelihood score to its physical characteristics.Our cost function measures the absolute average difference of the relative vorticity weighted by the distance between two consecutive time steps: where C is the cost function of a candidate track, N is the total number of the track's time steps, d is the distance between two consecutive track points and V is the relative vorticity at each time step.The number of potential tracks is quite large.However, their number can be significantly reduced by the application of a series of legitimate heuristics that remove tracks that present a "non-natural" behavior: (1) the location of the next candidate cyclone must be within a threshold range between successive time steps; (2) the maximum vorticity between the tracked cyclone and a candidate cyclone must not differ by more than 50 %; and (3) if the displacement between two successive displacements is more than 3 • long, then the angle between these displacements must be greater than 90 • .The first constraint prohibits the algorithm from searching for candidate features in the following time step in locations where the tracked cyclone could not be located.In our algorithm the cyclones are searched within a 5 • × 10 • latitudelongitude range.This is the largest possible displacement for extratropical cyclones as proposed by Hodges (1999).The second constraint prohibits the algorithm from choosing candidates that cannot be a possible evolution of the tracked feature.The use of a percentage is highly convenient since large vorticity values are subject to higher changes between consecutive time steps compared to smaller vorticity values.Finally, the third constraint prohibits the algorithm from taking into account a back-and-forth movement of the cyclone.Such displacements are more likely to take place in raw vorticity fields, where local maxima might change abruptly.For example, our algorithm would not choose the track passing from the points 2(12), 3(4) and 4(8) in Fig. 3, since the consecutive displacements present an angle of 74 • (marked in red), which is less than the 90 • threshold.
Finally, cycloTRACK returns one matrix as output for each track that contains information on the cyclone's track and its physical characteristics.Each matrix has one row for each of the track points and one column for each of the standard outputs and the optional physical diagnostics.These optional output diagnostics may vary depending on the study and the data inputs.Labeling the cyclonic circulations (Sect.2.1.1)and the cyclonic centers (Sect.2.1.2) within the tracks permits a post-processing analysis for determining potential merging and splitting of cyclones.For our application to the extratropical cyclones we consider only maximum 10 m wind speed and sea level pressure minima.As an example of the algorithm performance, Fig. 4 presents two cyclone tracks which evolve by sharing the same cyclonic circulation.Using the effective area diagnostic, in Fig. 5 we show the evolution of the two cyclones' relative vorticity, maximum 10 m wind speed and minima of sea level pressure.
It is likely that our method detects fronts associated with vorticity maxima as cyclone centers, especially when applied to high-resolution data sets, for example regional climatic simulations.In order to avoid the detection of a frontal zone, additional criteria of high or low complexity should be considered (e.g., Hewson and Titley, 2010).However, such criteria could be dependent on several factors -such as the spatial resolution of the data set -and would result in a "stricter" cyclone definition.The more precise the mathematical criteria, the more constrained the tracking results for systems of specific characteristics.For example, in the case of fronts, this could exclude the early stages of certain tracked cyclones that emerge from high-vorticity frontal areas of a "parent" cyclone.Figure 4 illustrates an example of a front detection.Inspection of surface pressure charts (not shown) showed that the first track point of the second cyclone (red dot in Fig. 4b) corresponds to the front of an extratropical cyclone (the one depicted by the black track).In the following time steps (Fig. 4c to f), this secondary vorticity maximum evolves into a strong cyclone (red track) which presents its own low-pressure minimum.Here we capture the initial stage of the vorticity maximum, before the occurrence of a pressure minimum.Nevertheless, skipping the application of additional criteria may require a post-processing of the resulting tracks in order to exclude the "wrong" ones or those that do not match the research needs.

Application of cycloTRACK in a climatological context and analysis of its sensitivity
In this section we present the results of the application of our algorithm for all winters (December-January-February) of the period 1989-2009.We also present the results on three sets of sensitivity tests: (a) on relative vorticity filtering, (b) on the cost function of Eq. ( 2), and (c) on the constraint that relative vorticity between two consecutive track points must not differ by more than 50 %.In all sensitivity tests, the threshold used to define cyclones is 3 × 10 −5 s −1 , and we analyze only those tracks with lifetimes longer than 1 day.

Sensitivity of filtering the relative vorticity field
We tested three different filter strengths (described in Sect.2.1.1)with the ERA-I data set.The applied spatial filters correspond to 3 × 3, 5 × 5 and 7 × 7 grid point filtering, named filter3, filter5 and filter7, respectively.Figure 6a presents the number of detected cyclonic centers as a function of their relative vorticity for all three sensitivity tests, and Fig. 6b presents their relative frequency.Since all tests are bounded to identify cyclones exceeding a common threshold of 3 × 10 −5 s −1 , and since filtering decreases the relative vorticity values due to its smoothing operation, it is not surprising that the total number of detected cyclone centers is reduced with increasing filtering intensity.Regardless of the spatial filtering strength, all three sensitivity tests present an exponential distribution (Fig. 6a), and the stronger the filter, the more the cyclones' intensities are reduced (Fig. 6b).
Strong filtering versus weak filtering may have two effects.First, it tends to detect fewer tracks which also correspond to the stronger cyclones.Second, it tends to reduce the cyclone track lengths by not taking into account the weakest vorticity perturbations in the early and late stages of a cyclone track.The validity of the first hypothesis is evident from Fig. 1, where smoothing suppresses many weak cyclonic centers, but stronger cyclones (such as Anatol) are equally detected with all three filters.To verify the second hypothesis, we investigate the characteristics of the tracks as detected by fil-ter3, filter5 and filter7.Figure 7a, b and c show the distribution of the relative frequency for the lifetime of cyclone tracks, the average speed of the cyclones and their maximum relative vorticity.We observed no significant changes on the results obtained with the different filters when considering the cyclone lifetime.Consequently, the second hypothesis, that average track characteristics are sensitive to filtering, can be rejected.It is interesting, however, that our applications using weak filtering detect weak cyclones that have similar lifetimes.The distributions of the relative frequencies of average speed of cyclones are similar for all three filters (Fig. 7b).This means that the weaker cyclones in filter3 and filter5 do not correspond to weak stationary vorticity perturbations, and they also do not evolve into strong extratropical cyclones.The reasons for not evolving into strong cyclones is an interesting issue; however investigation of this is beyond the scope of this paper.Figure 8 shows the cyclone center density (CCD) for all three filtering strengths.It is evident that different magnitudes of CCD are observed that depend on the filtering strength; however, the spatial pattern remains coherent for all three cases.A question that may arise is whether all sensitivity tests share the same cyclone centers, while additional weak centers are detected in filter3 and filter5 that are suppressed in filter7 due to the smoothing operation.To address this question we took into account all points of the distributions in Fig. 6 and we associated the common points between filter3 and filter7 with one another (points sharing the same timing and having a distance inferior to 5 • ).Results showed that filter7 shared 52 % of its points (2331 points) with filter3.The median of the intensity of the common points of filter3 corresponded to the 78th percentile of all filter3 points' intensity.Consequently, cyclones in filter7 correspond to the strongest cyclones of the weakly filtered data.This is in accordance with the relative frequency of cyclone centers' intensity in Fig. 6b, where most identified cyclones using filter7 are concentrated as weaker relative vorticity values compared to filter3 and filter5.
The effect of filtering, for example filter7 compared to fil-ter3, is characteristic to the CCD within the Mediterranean region, where the cyclones are known to be weaker than the extratropical cyclones forming over the oceans ( Čampa and Wernli, 2012).Indeed, in filter7 there is a dramatic decrease in detected cyclones over the Mediterranean Sea compared to filter3 and filter5.Figure 8 presents a high similarity with the results of other algorithms (Neu et al., 2013), independent of whether a filtering was applied or whether the sea level pressure and the relative vorticity were used as inputs for the detection of cyclones.Indeed, CCD maxima are distinctly located over the Pacific Ocean, the North Atlantic Ocean and the Mediterranean.Furthermore, regardless of the filtering strength, the relative frequency distributions of cyclones speed and lifetime (Fig. 7a and b) seem to be in good agreement with the other algorithms (Neu et al., 2013) presenting most probable cyclone speeds between 30 and 40 km h −1 and cyclone lifetime relative frequency distributions decreasing  exponentially from less than 2 days up to a total of approximately 8 days.
Figure 9 presents the annual number of cyclone centers.For all three filters, our results are in agreement with those of Neu et al. (2013), showing no specific interannual trend.As expected, the number of cyclone centers per year depends on the filtering strength.It decreases from approximately 9000 per year for filter3 to approximately 3000 per year for filter7.All three tests are within the ranges of other algorithms that range from 2000 to 12 000 per year.However, it is only fil-ter5 that is consistent with the majority of the results of other algorithms, which calculated 4000 to 7000 cyclonic centers per year.The time series phasings are in good agreement between filter3 and filter5, presenting a correlation of 0.91.On the other hand, correlation between filter5 and filter7 is 0.43, suggesting that the time series phasing between the two sensitivity tests is dependent on the weaker cyclones that are suppressed in filter7.This should not raise a question as to the soundness of the different test results but rather as to the sensitivity of the results to the different filtering strengths.

Sensitivity of tracking parameters
We performed two additional sets of sensitivity tests in order to assess the efficiency of our tracking method (step 2).The first set relates to the cost function (Eq.2) and it is composed of (a) S rel , where the final track choice in step 2 is only dependent to the track relative vorticity evolution (Eq.3), and (b) S dist , where the cost function is only dependent on the distance between consecutive track points (Eq.4).The second set of sensitivity experiments relates to the constraint that the relative vorticity between consecutive track points may not vary by more than 50 % (Sect.2.2) and it is composed of three tests, where the 50 % threshold has been modified to 25 % (S 25 % ), 75 % (S 75 % ) and 100 % (S 100 % ).In all three tests we used the original cost function (Eq.2) and the identified cyclones from filter3.We focused on filter3 since this is the data set with the highest number of identified cyclones (Fig. 6), thus amplifying the differences between the tracking results of the sensitivity tests.
Figure 10 presents the lifetime of tracks (i.e., the number of track points) and their average speed (i.e., distance between the track points) for both sets of sensitivity experiments.The results of the first set of experiments that focus on the cost function (Fig. 10a and b) show that the cyclones' lifetime and average speed is similar for all sensitivity tests: filter3, Srel and Sdist (maximum differences are less than 1 %).This suggests that the lifetime of track points and their average speed are rather insensitive to the change of the cost function.This is due to the fact that the algorithm always presented several alternative tracks for a single cyclone but that, in the majority of the cases, these alternative tracks were similar and only presented small deviations from the cyclones' main path.In such cases, the usefulness of the cost function is in choosing the smoothest track in terms of intensity and distance between consecutive track points.It is worth noting that, in S rel and S dist , the algorithm was still bounded by the constraint of linking cyclone centers that presented relative vorticity values that did not vary by more than 50 %.Climatologically, the term d in the cost function does not significantly change the results of the algorithm.However, this term was shown to play an important role, especially when several cyclones of quasi-equal vorticity were located within the 10 • × 5 • area but unrealistically far from the tracked center.The results of the second set of experiments that relate to the 50 % threshold (Fig. 10c and d) reveal similar distributions for all varying thresholds; however when comparing S 100 % and S 25 % , the former tends to form longer tracks (Fig. 10c) with longer distances between the track points (Fig. 10d).Indeed, when applying smaller (higher) thresholds to the permitted evolution of the cyclones intensity, it is more likely that the algorithm will form shorter (longer) tracks due to the smaller (higher) accepted differences on the relative vorticity evolution of consecutive track points.Ideally, the 50 % threshold could be neglected.However, this would create numerous alternative tracks in the cases of highresolution data sets.In general, the constraints applied in step 2 (i.e., 50 % threshold, searching cyclones within a 10 • × 5 • area and the angle criterion; Sect.2.2) have been found to be a fair compromise between cutting off "unnatural" candidate cyclone tracks and providing all possible tracks for the algorithm to depict the "correct" one according to the cost function.

Physical coherence of the tracked cyclones
In this section we perform an analysis of the effective area diagnostic tool described in Sect.2.1.3,taking into account only filter5 results.Figure 11 presents the composite life cycle of the cyclones physical characteristics, centered on the time of the maximum vorticity of the tracks (mature stage) and averaged for all tracks detected in the Pacific Ocean (from 130 to 240 • longitude and from 30 to 90 • latitude), North Atlantic Ocean (from 300 to 360 • longitude and from 30 to 90 • latitude) and within the Mediterranean region (from 345 to 45 • longitude and from 25 to 50 • latitude).The results show that, regardless of the region, there is a strong coherence between the life cycle of sea level pressure minima, relative vorticity and maximum 10 m wind speed.Cyclones tend to intensify rapidly but weaken at a slower rate.For the construction of the composites, there is no distinction regarding the cyclones' lifetime.Also, we should note that the further we get from the time of the cyclone's maximum vorticity (i.e., the composite center), the fewer cyclones that last long enough to provide diagnostics for the composites.For example, the Mediterranean cyclones' lifetime is inferior to the other extratropical cyclones and rarely exceeds 2-3 days.Nevertheless, our motivation here is to assess the validity of the effective area diagnostic, which seems to correctly capture the physical characteristics of the life cycle of cyclones regardless of the region.Indeed, in agreement with Čampa and Wernli (2012), we found that Mediterranean cyclones are less deep, in terms of sea level pressure, while Atlantic cyclones are slightly deeper than those occurring over the Pacific Ocean.

Conclusions
In this article we presented a new algorithm for identifying and tracking cyclones, applied to winter extratropical cyclonic systems over the Northern Hemisphere.We tested our algorithm's performance for different filtering strengths applied to the high-frequency relative vorticity fields.Results showed that the number of tracks was inversely proportional to the filter strength, while the cyclone spatial and temporal variability were coherent with those produced by other tracking algorithms presented in the literature.Finally, our algorithm successfully captured the physical characteristics of cyclones.
Our identification and tracking algorithm uses as few constraints as possible, not only for tracking weak vorticity perturbations which evolved into strong cyclones but also for tracking weak perturbations which did not evolve into strong cyclones.This allows for not only assessment of the algorithm's sensitivity to data filtering but also, in the future, a more precise description of the environmental conditions which favor cyclone intensification.Furthermore, we chose the vorticity criteria to vary dynamically (vorticity must not vary by more than 50 % in consecutive time steps) and we avoided any threshold or cut-off values.It should be noted that although we applied our algorithm based on relative vorticity to identify and track cyclones, the same algorithm might be applied to any data set which presents enclosed areas after applying a threshold value.For instance, our algorithm could be applied for tracking supercells or mesoscale convective systems using data sets of brightness temperature or cloud cover.
Our tracking approach is based on minimizing a cost function of vorticity maxima.We observed some mistakes, especially when cyclonic circulations were found to be very noisy with multiple local maxima.As an alternative cost function, it would be interesting to explore the weighted mean differences of additional cyclone physical characteristics (pressure, wind speed etc.) between consecutive time steps.This has been previously applied by Machado et al. (1998) for tracking MCSs based on brightness temperature satellite observations.However, their method assumes a priori choice of the weighting value, risking restraining our method's adaptability in tracking cyclones of different origin (e.g., extratropical and tropical cyclones).Our algorithm links cyclone centers in consecutive time steps, in contrast with the alternative configuration proposed by Machado et al. (1998) andInatsu (2009) to link enclosed areas.This decision was made because large cyclonic circulations would not correspond to a single cyclone if enclosed areas were linked, and additional criteria -and/or filtering -would be needed, while weak cyclones would be neglected.
Further development of our algorithm includes the extension of the identification part to three dimensions and the extension of the method's adaptability for different atmospheric features such as MCSs.CycloTRACK was implemented in MatLab, and its source code is freely available from the corresponding author upon request.

Figure 1 .
Figure 1.(a) Relative vorticity raw fields at 00:00 UTC, 3 December 1999.The applied threshold is 3 × 10 −5 s −1 .Crosses represent central maxima located in the center of a 3 × 3 grid point area.(b) As in (a) but relative vorticity field is filtered using a 3 × 3 correlation spatial filter.(c) As in (a) but relative vorticity field is filtered using a 5 × 5 correlation spatial filter.(d) As in (a) but relative vorticity field is filtered using a 7 × 7 correlation spatial filter.

Figure 2 .
Figure 2. The storm "Anatol" at 00:00 UTC, 3 December 1999.Relative vorticity is smoothed by a 3 × 3 spatial filtering (color bar), the contour denotes the mean sea level pressure, and arrows denote 10 m wind field.The thick black circle represents the effective area of the cyclone.Locations and values of maximum wind speed and lower sea level pressure are depicted by the thick lines.

Figure 3 .
Figure 3. (a) An idealized case of cyclone locations in four time steps.The locations are depicted by circles.The numbers above the locations are in the form X(Y ), where X denotes the time step and Y denotes the relative vorticity.The circle sizes are proportional to the relative vorticity value of the cyclone.(b) All possible trajectories of cyclone 2(12) by searching backward and forward in time.(c) A selected subset of tracks of (b): we retain those which present the minimum average change of relative vorticity in successive time steps.

Figure 4 .
Figure 4. Relative vorticity smoothed by a 3 × 3 spatial filter (color bar); sea level pressure (contours, with a 5 hPa interval; thick contour denotes 1000 hPa); and tracks (thin lines) of two splitting cyclones for different times during December 1999.

Figure 5 .
Figure 5. (a) Maximum relative vorticity (solid line) at the track centers and minimum sea level pressure (dashed line) as detected within the cyclones effective area for the two cyclones shown in Fig. 4. (b) As in (a) but dashed line corresponds to maximum 10 m wind speed.Time series coloring is consistent with the tracks colors in Fig. 4. The horizontal axes represent the period 6-16 December 1999.

Figure 6 .
Figure 6.(a) The number of cyclonic centers as function of their relative vorticity and as detected by the three algorithm sensitivity tests.(b) Relative frequency distributions as a function of the relative vorticity for the identified cyclone centers.

Figure 7 .
Figure 7. (a) Relative frequency distributions of cyclones lifetimes for the three sensitivity tests.(b) As in (a) but for cyclones average speed.(c) As in (a) but for tracks' maximum relative vorticity.

Figure 8 .
Figure 8. Cyclone center density expressed as the percentage of cyclone occurrence per time step and per unit area of 1000 km 2 for (a) filter3, (b) filter5 and (c) filter7.

Figure 9 .
Figure 9.The number of cyclone centers per year for the three sensitivity tests.

Figure 10 .
Figure 10.(a) Relative frequency distribution of cyclones' lifetime for the sensitivity tests filter3, S rel and S dist .(b) As in (a) but for cyclones average speed.(c) As in (a) but for the sensitivity tests filter3, S 25 % , S 75 % and S 100 % .(d) As in (b) but for the sensitivity tests filter3, S 25 % , S 75 % and S 100 % .

Figure 11 .
Figure 11.(a) Average composite time series of Pacific cyclones physical characteristics."0 h" corresponds to the time that the cyclone presents its maximum relative vorticity.We denote relative vorticity (thick black line), sea level pressure (red thick line) and maximum 10 m wind speed (thin black line).Wind speed scale values are shown in the left vertical axes in parentheses.(b) As in (a) but for the Atlantic cyclones.(c) As in (a) but for Mediterranean cyclones.Note that y axis does not have the same value intervals across the three panels.