{'fraction_non_alphanumeric': 0.026592586309271353, 'fraction_numerical': 0.015056408516413607, 'mean_word_length': 4.923708232141362, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A 4D feature-tracking algorithm: A multidimensional view of cyclone systems

4500. An objective four-dimensional (4D) algorithm developed to track extratropical relative vorticity anomaly 3D structure over time is introduced and validated. The STACKER algorithm, structured with the TRACKER single-level tracking algorithm as source of the single-level raw tracks, objectively combines tracks from various levels to determine the 3D structure of the cyclone (or anticyclone) events throughout their life cycle. STACKER works progressively, beginning with two initial levels and then adding additional levels to the stack in a bottom-up and/or top-down approach. This allows an iterative stacking approach, adding one level at a time, resulting in an optimized 4D determination of relative vorticity anomaly events. A two-stage validation process is carried out with the ECMWF reanalysis ERA-Interim dataset for the 2015 austral winter. First the overall tracking capability during an austral winter, taking into account a set of climate indicators and their impacts on Southern Hemisphere circulation, was compared to previous climatologies, in order to verify the density and distribution of the cyclone events detected by STACKER. Results show the cyclone density distribution is in very good agreement with previous climatologies, after taking into account potential differences due to climate variability and different tracking methodologies. The second stage focuses on three different long-lived events over the Southern Hemisphere during the winter of 2015, spanning seven different pressure levels. Both GOES satellite imagery, infrared and water vapour channels, and ERA-Interim cloud cover products are used in order to validate the tracks obtained as well as the algorithm’s capability and reli-ability. The observed 3D cyclone structures and their time evolution are consistent with current understanding of cyclone system development. Thus, the two-stage validation confirms that the algorithm is suitable to track multilevel events, and can follow and analyse their 3D life cycle and develop full 3D climatologies and climate variability

An objective four-dimensional (4D) algorithm developed to track extratropical relative vorticity anomaly 3D structure over time is introduced and validated. The STACKER algorithm, structured with the TRACKER single-level tracking algorithm as source of the single-level raw tracks, objectively combines tracks from various levels to determine the 3D structure of the cyclone (or anticyclone) events throughout their life cycle. STACKER works progressively, beginning with two initial levels and then adding additional levels to the stack in a bottom-up and/or top-down approach. This allows an iterative stacking approach, adding one level at a time, resulting in an optimized 4D determination of relative vorticity anomaly events. A two-stage validation process is carried out with the ECMWF reanalysis ERA-Interim dataset for the 2015 austral winter. First the overall tracking capability during an austral winter, taking into account a set of climate indicators and their impacts on Southern Hemisphere circulation, was compared to previous climatologies, in order to verify the density and distribution of the cyclone events detected by STACKER. Results show the cyclone density distribution is in very good agreement with previous climatologies, after taking into account potential differences due to climate variability and different tracking methodologies. The second stage focuses on three different long-lived events over the Southern Hemisphere during the winter of 2015, spanning seven different pressure levels. Both GOES satellite imagery, infrared and water vapour channels, and ERA-Interim cloud cover products are used in order to validate the tracks obtained as well as the algorithm's capability and reliability. The observed 3D cyclone structures and their time evolution are consistent with current understanding of cyclone system development. Thus, the two-stage validation confirms that the algorithm is suitable to track multilevel events, and can follow and analyse their 3D life cycle and develop full 3D climatologies and climate variability studies.
While tropical cyclones, often referred to as hurricanes or typhoons depending on their location, have strong seasonality and characteristic tracks in the Tropics, extratropical cyclones have different life cycles and displacements, are much more frequent and can vary both in shape and structure spanning a diverse range of sizes. Ayrault and Joly (2000) identified 12 cyclone families over the extratropical North Atlantic in a study of approximately 6,000 extratropical cyclone events. Furthermore, extratropical cyclones can develop in diverse synoptic situations at mid-to high latitudes depending on their relationship with the evolution of the jet stream, polar ice coverage and midlatitude sea-surface temperature (SST) gradients. Some cyclones are confined to lower tropospheric levels and others may develop to extend through the troposphere (e.g. shallow and deep cyclones: Lim and Simmonds, 2007). Chen et al. (2014) found that the cyclone and anticyclone activity over China are linked to upper tropospheric jets over East Asia.
The detection of cyclone tracks is a challenge that dates back to the nineteenth century: its origin can be found in van Bebber (1891). Petterssen (1956) studied the occurrence and behaviour of moving low-pressure centres in the Northern Hemisphere. Subsequent analyses were carried out by several authors in the same region, focusing on climatologies of cyclone development (Homar et al., 2003) and their relationship with atmospheric dynamics (Hoskins and Valdes, 1990). To efficiently study and to quantify cyclone systems, a number of algorithms have been proposed in recent years, especially considering the renewed interest in the context of climate change. The Intergovernmental Panel on Climate Change has noted the relevance of understanding the dynamic climatology of cyclones and its relation to climate change in several assessments (IPCC, 2007;2014). Grieger et al. (2018) also highlight the relevance of such studies over the Southern Hemisphere (SH) given the crucial role extra cyclones play in the poleward transport of heat and moisture to polar regions.
One of the most frequent ways to assess the cyclonic activity involves statistical criteria, which in general are focused on the occurrence and position of the systems, but not generally on the tracking or their physical properties (e.g. Jansa et al., 2001;Finnis et al., 2007). A more specific direct approach to identify and analyse these systems is to use objective feature-tracking algorithms that detect features belonging to cyclones in gridded reanalysis products. A number of cyclone tracking algorithms have been developed using different features to detect and track cyclones and anticyclones (Neu et al., 2013 and references therein;Grieger et al., 2018 and references therein; Hodges, 1995;1999;2008). A widely used tracking method, the neighbour point tracking (NPT), consists of a two-step approach in which first, the cyclonic system centres, defined as a local maximum or minimum value of a climatic variable, are detected at each time step and then the track based on the nearest-neighbour concept (Kelemen et al., 2015) is obtained through the temporal connection between cyclone centres (Murray and Simmonds, 1991a,b;Satake et al., 2013;Flaounas et al., 2014). Most frequently used climate variables for tracking purposes include mean-sea-level pressure (MSLP), MSLP gradients, geopotential height, potential and relative vorticity, singly or in combination (hybrid approach), e.g. Trigo et al., 1999;Jansa et al., 2001;Picornell et al., 2001;Hoskins and Hodges, 2002;Hanson et al., 2004;Bartholy et al., 2008;Hewson and Titley, 2010. Such algorithms imply the use of a decision tree in order to apply necessary constraints regarding the detection and evolution of the feature under study. Thus, a robust and efficient tracking algorithm has to be able to identify feature values, usually below or above a defined threshold, their path over time, and also the link between the detected centres as well as their life cycle. But in addition, the algorithm also has to consider that cyclones can split or merge with other cyclones. Grieger et al. (2018) carried out an assessment on the degree of consensus among 14 state-of-the-art cyclone identification and tracking methods, following Neu et al. (2013). Both studies found that the selected feature-tracking algorithms yielded similar overall climatic structures, but notable differences remain, e.g. in intensity of events and event counts. Tracking results tend to group into vorticity-based and MSLP/geopotential-based clusters. Different cyclone tracking schemes yielded different cyclone trends, sometimes differing even in sign.
Currently available algorithms mostly track atmospheric events on a single geopotential/pressure level, a few of them considering vertically averaged levels, multivariable ensemble or models (Hodges et al., 2003;Serra et al., 2010;Kelemen et al., 2015). Because cyclones can have diverse characteristics and evolution at various levels, single-level or average tracking may impact climate analysis. For example, over the SH the Andes mountain range and the African Plateau in the SH (Inatsu and Hoskins, 2004;Hoskins and Hodges, 2005), different choices of tracking level can give different results for storm track and cyclone evolution. Thus determining cyclone track structure is an important issue for climate analysis. Lim and Simmonds (2007) compared cyclone tracks at six levels, between MSLP and 500 hPa, to assess the vertical climatology of SH cyclones for the austral winter months June-August (JJA), for the period 1979-2001. Detected cyclone positions at the various levels were combined by considering the distance between cyclone centres at a given time, as given by the pressure anomaly (depth) at MSLP and geopotential height anomaly centres, from one level to the next above, using a reference radius of 444 km. Starting at sea level, cyclone centres detected at two adjacent levels, within this reference radius, were considered as belonging to the same cyclone event. At each time step this vertical coupling was carried out up to the level that complied with this criterion. Cyclone tracks were then defined in terms of multiple-level events, by considering their variations in height, at the times of overlap only. Results of this bottom-to-top approach composition were further separated into shallow (up to 700 hPa) and well-organized cyclones (at least up to 500 hPa) in an SH cyclone climatology study. Extending and improving the track coupling process to several more pressure levels would provide new insights on cyclone types and cyclone life cycles as well as help assess linkages with other weather and climate processes. In addition, a multilevel tracking developed in this way should take into consideration the different track lengths at various levels, including track points existing at different levels, prior to or after a period of vertical overlapping. In other words, the description of the multilevel system should not be restricted only to the period of strict overlapping at different heights, in order to fully describe the system's evolution.
The method developed by Lim and Simmonds (2007) provides the starting point for the implementation of a more comprehensive approach to the processing and combination of cyclone tracks at multiple levels for detailed studies. First, it is necessary to extend the detection process to a larger number of pressure levels throughout the troposphere. Secondly, given the different kinds of cyclone and anticyclone processes and life cycles, it is necessary to apply both a bottom-to-top approach as well as a top-to-bottom one. Thirdly, it is necessary to extend the climatology to span recent years as well as all seasons. This implies the processing of vast amounts of data for the detection of tracks at each level and the subsequent combination of numerous tracks from levels, both bottom-to-top and top-to-bottom, to determine the vertical structure and space-time evolution of events. This requires the implementation of algorithms that use an optimized, automated and efficient approach to assess and define the vertical structure of cyclones and anticyclones. From now on, reference will exclusively be made to cyclones, though the method may also be used for anticyclones.
The present study is based on the cyclone tracking code developed by Hodges (1995;1999) which uses relative vorticity rather than pressure or geopotential anomalies computed from MSLP and geopotential as in Keable et al. (2002) or Lim and Simmonds (2007), the latter being based on geostrophic vorticity computed from the geopotential. Following a dynamic programming approach, an automated algorithm, called STACKER, was developed to determine which cyclone events detected at the different levels are to be combined into single events, preserving all the track features detected at each level. The STACKER algorithm is built using a smoothness cost-function approach. European Centre for Medium-range Weather Forecasts (ECMWF) ReAnalysis ERA-Interim products for the SH were used. Not only are there fewer studies for the SH, but the oceanic nature of the hemisphere and major orographic barriers -such as the Andes mountain range, the well-established regions of cyclogenesis on the eastern seaboard of southern South America and cyclolysis on its western seaboard, on either side of the Andes mountain range (Inatsu and Hoskins, 2004;Hoskins and Hodges, 2005) -provide interesting tests for the algorithm and its validation.
The STACKER algorithm is introduced in section 2. Multilevel cyclone tracking for the 2015 winter (JJA) was carried out with STACKER. Section 3 presents the validation of the STACKER algorithm using the 2015 JJA mean state and cyclone density. Furthermore, three selected multilevel events developing near southern South America are introduced and compared with GOES imagery, ERA-Interim cloud cover and winds for further algorithm validation. Conclusions are presented in section 4.

Data
In this section the 4D STACKER algorithm, developed to track significant extratropical anomaly structures by using relative vorticity anomalies and multilevel track coupling, is presented. The reanalysis data product used to help develop and test the algorithm is the ERA-Interim 6-hourly relative vorticity (RV), with a horizontal resolution of 1.5 • × 1.5 • at various levels between the surface and 100 hPa (lower stratosphere). In terms of choice of ERA-Interim for the testing of STACKER, it should be noted that Hodges et al. (2003) compared four different reanalyses (Japanese JR25, ERA-Interim, National Aeronautics and Space Administration Modern Era Retrospective -analysis for Research and Applications (NASA MERRA) reanalysis and the National Centers for Environmental Prediction Climate Forecast System Reanalysis (NCEP CFSR)) for composite cyclone diagnostics and showed that the detected cyclones compared well between the reanalyses in terms of numbers, location and spatial distribution, with differences being primarily found in the cyclone intensities. Froude (2010), working with a set of reanalysis resolutions, concluded that the relative performances of the different Ensemble Prediction Systems (EPS) in predicting cyclones seem to agree with the results obtained using more conventional forecast verification measures (Park et al., 2008), and results comparing the skill of the ensemble mean with that of the control forecasts show that for cyclone position there is very little difference between all of the EPSs, independent of resolution. All tests with STACKER were carried out over the Southern Hemisphere during the winter season (JJA) of 2015.

Tracking algorithm
The base algorithm selected to track the cyclones at different levels in the present study is the TRACK algorithm (Hodges, 1995) the details of which can be found in several papers (Hodges, 1994;1995;1999;2008 Dacre et al., 2012;Zappa et al., 2013;2014;Pinheiro et al., 2017) and its implementation is well documented. Grieger et al. (2018) also conclude that vorticity-based tracking algorithms can identify cyclone structures at an early stage in their life cycle, hence their use in the current implementation of the STACKER algorithm. A brief review of TRACK is included here in order to understand the subsequent steps. TRACK was designed to track feature points in a gridded dataset over consecutive time steps on a single level. The first step of the algorithm spatially filters the data to remove the planetary scales and to reduce the noise in the RV field to allow more reliable tracking of the features (Hoskins and Hodges, 2002). The filter band is chosen for the appropriate scales of the feature of interest. This is done by transforming the field to spectral space using spherical harmonics, applying a triangular truncation filter and transforming back to real space. For synoptic-scale cyclones the filter band chosen is T5-42, which has been used in the studies referenced above. A tapering filter is also applied to spectral coefficients to reduce Gibbs oscillations (Hoskins and Sardeshmukh, 1984). The efficacy of this approach is highlighted in Froude (2010) (their figure 1) which compares the T42 RV filtered field and the feature points identified with the much higher resolution T799 (25 km) MSLP field. When using RV for synoptic-scale cyclones the same filtering is performed irrespective of the original higher data resolution.
The next step of the algorithm identifies local minimum or maximum values of RV as feature points representing cyclonic and anticyclonic events in the SH. These feature points are then formed into tracks by first initialising a set of tracks based on a nearest-neighbour method and then using the minimization of a Smoothness Cost Function (SCF) subject to adaptive constraints (Hodges, 1999) suitable for extratropical cyclones, to produce an optimal smoothest set of tracks. This whole procedure can be summarized in three main steps: Data pre-processing, Feature Point Detection, and Tracking.
Previous studies have compared the automated tracking using TRACK with a more subjective manual tracking approach, though this can be sensitive to different observers. Colle et al. (2013) found a success rate ∼85% using RV over a North American-western Atlantic region that improves if short-lived systems, ∼1 day, are excluded. The success rate was found to increase to ∼90-95% using MSLP. The study by Gramcianinov, Hodges and Camargo ("The properties and genesis environments of South Atlantic cyclones," submitted, 2018), using RV, found a success rate over a region of the South Atlantic of ∼89%. These success rates compare favourably with other schemes.

STACKER algorithm
The STACKER algorithm is designed so as to associate all tracks obtained from a given feature-tracking algorithm that represent the same synoptic event at different vertical levels; here pressure levels are used. In order to expand the analysis to include more pressure levels and hence to describe the four-dimensional evolution of cyclones a binary linear programming model is developed, in which each possible combination of simultaneous detections on adjacent pressure levels would become an optimization variable to represent the most coherent relationship between the given set of detections. This algorithm named "STACKER algorithm" is a progressive model optimization algorithm. It considers tracks from TRACK output at two adjacent pressure levels and links them according to an objective function criterion. Subsequently, the resulting stacked tracks are linked with the next adjacent pressure-level tracks, until the stacking is complete. The name STACKER comes from the idea of stacking tracks from upper levels over tracks from lower levels. The algorithm works in either direction, bottom-to-top and top-to-bottom, i.e. including the possibility of stacking lower-level tracks under higher-level tracks, such as upper-level troughs or cut-off lows. STACKER only works with track detections from two levels at a time and iterates over the levels adding a single level to the stack each time. Figure 1 presents the STACKER algorithm flowchart. The stacking process is highlighted in blue, while the single-level tracking input is highlighted in green. The flowchart shows the decision-making process and the optimization process, given all the possible outcomes of the decision process. The optimization process should create a link for every possible pair of combinations between two consecutive levels, i.e. track-to-track linking. However, this approach assumes that each lower track corresponds to only one upper track, but this may not be entirely true. Indeed, it is possible that an event changes during the evolution of its life cycle as a consequence of changes in atmospheric circulation, land-sea transitions, or it may be impacted by significant orography, resulting either in lysis or, occasionally, modification of part of the vertical structure while maintaining its nature. To avoid this limitation a point-to-point linking is used. This method considers a link for every possible combination of a lower (upper) point with an upper (lower) point, while constraining the problem such that a link between two points can be made only if all of the other links can be made simultaneously. Furthermore, points must not be associated with more than one active link. Thus the optimization process returns a track-to-track association born from a point-to-point correspondence subject to chosen constraints ( Figure 1). Specific details of the STACKER algorithm can be found in the appendix.
Using this approach, the algorithm maintains the coherence of linking two potentially associated tracks assuming they belong to the same object and thereby gains the flexibility of linking more than one track on each pressure level. Nonetheless, the combinatorial problem becomes far more complex and the amount of optimization variables, problem constraints and processing resources needed increases considerably. In order to avoid highly improbable associations between very dissimilar detections, a "linkage" threshold value, i.e. a strong constraint, will discern which associations will be relevant for the study and thereby decrease the problem complexity. Also, some filters can be applied, like maximum radius distance between points, minimum time duration, minimum distance travelled, region of origin, etc. Figure 2 shows how a maximum radius threshold is considered at each point along a track to constrain the coupling of points from the upper track (blue disk) in a manner similar to that used by Lim and Simmonds (2007). The plot also highlights how this helps discriminate the coupling process for tracks crossing each other in the lower level.
Due to the algorithm's stacking nature, which operates in both directions, bottom-to-top and top-to-bottom, from now on and for the sake of simplicity the section that has already been stacked or is the result of a previous stage of the algorithm will usually be referred to as "lower", independently of the direction of the stacking process. Similarly, what it is commonly referred to as "upper" is the new section that is to be added or stacked onto an already existing stack. Figure 3 explains the iterative stacking process through a simple example considering track segments detected at three adjacent levels. "Upper" tracks are usually one-level tracks, while "lower" tracks are considered to contain tracks from more than one pressure level, resulting from the stacking process, i.e. stacks, except during the initial stack iteration where one "lower" and one "upper" track are combined. While the same process can be used in either stacking direction, bottom-to-top or top-to-bottom oriented, the final solution, while similar may not exactly be identical for a given multi-levelled event stacked in different directions. Minor differences may result from potential fragmentation of tracks ( Figure 3). The differences between both directions could be attributed to the physical nature of the cyclonic system. Thus the best direction would be associated with what physical processes are driving cyclone development or maintenance. For instance, at upper levels in subtropical latitudes, cyclogenesis could be due to baroclinic instability, while at lower levels at polar latitudes, sea-ice to open water thermal contrast could contribute to the forcing at lower levels (Hakim, 2003;Joly et al., 2003 and references therein). Bottom-to-top and top-to-bottom tracks which overlap both in space and time are merged to obtain an integrated view of the event. This merger feature has been included here as an automatic process given the use of STACKER for climatology studies (see further details in appendix). Studies focusing on case-study dynamics may consider bottom-to-top and top-to-bottom, presenting results independently if desired.
Before proceeding further, it is important to note an issue linked to the relative vorticity threshold value used in track determination. While this threshold is defined for TRACK, it is also relevant to STACKER, using a vorticity-based feature-tracking, because its choice may impact the stacking process. This happens because tracks obtained at various levels are combined in the stacking process. The relative vorticity threshold as defined in TRACK is fixed at −1.0 × 10 -5 s −1 , without height or latitude dependence. This is set at a comparatively low value to capture as much of the cyclone life cycle as possible for any season. A more restrictive threshold will reduce the number of systems identified but will also shorten the life cycles of those that are identified, more so for the weaker summer systems. Furthermore, it is important to note that TRACK, as well as other vorticity-based tracking algorithms (Grieger et al., 2018), spectrally filters the vorticity to remove the small-scale noise (see section 2). This will reduce the tendency to identify a large number of small-scale systems at all levels. RV exhibits a somewhat larger mean value at 500 hPa, the non-divergence level. Inspection of monthly mean RV fields shows that the range of RV values does not change significantly with height, nor with season. RV threshold levels for the lowermost and uppermost levels are similar. Furthermore, variations are primarily spatial. RV standard deviation fields show that the amplitude of RV perturbations increase with height. Hence a fixed threshold at a low value does not significantly impact the detection of cyclonic systems.
Furthermore, a fixed RV threshold value such as the one chosen here does not affect cyclone event counts in the lower levels near the surface nor in levels near the tropopause. When stacks are constructed by STACKER bottom-up or top-down and the stacks contain tracks near the upper or lower levels this does not pose a problem because those initial levels, with similar threshold values, will constrain the number of tracks from the middle troposphere which are stacked. Concern about this value could arise for stacks which begin at or around 500 hPa, i.e. mid-tropospheric systems, where this threshold may enhance the number of tracks detected there, since in principle relative vorticity maximizes there. If this is so, and if the threshold value becomes more restrictive, the number of systems identified would be reduced, and, as already noted, result in shorter life cycles for those identified. This criterion has precedents in, for example, Pinheiro et al. (2017). Even so, and in order to test the threshold sensitivity, the TRACK algorithm was run at 500 hPa, with a set of more restrictive RV threshold values. Results indicated that the difference may be limited to less than 2% fewer detected events, which is negligible, even more so if the discrepancies between the different tracking algorithms compared in Neu et al. (2013) and Grieger et al. (2018) are taken into consideration. Note that a similar issue could arise with other tracking algorithms when used with STACKER, depending on the criteria used to define the anomalies to be tracked at different levels, even with pressure or geopotential-based tracking.
The detailed definitions and concepts of essential aspects, as well as the mathematical components of the STACKER algorithm, are available in the appendix.
The STACKER algorithm only uses detected cyclone tracks as input, i.e. cyclone position at a given time and height, as the above paragraphs show. The issues regarding RV pertain to the use of the RV-based TRACK algorithm. STACKER, however, can be used in principle with tracks obtained by any reliable cyclone tracking algorithm. As was the case here, given the RV threshold dependence with height, care will be required to determine the adequate height-dependent threshold for the variable used in the selected tracking algorithm, be it pressure, geopotential, temperature or a combination thereof.

VALIDATION STUDIES
The fields used for testing the combined tracking and stacking algorithms came from the ERA-Interim reanalysis dataset. As explained above, each event was first independently tracked with TRACK on each pressure level during the 2015 JJA (austral winter) period and then stacked with the STACKER Algorithm (SA). A total of 753 multilevel, stacked extratropical cyclone events over the SH (south of 30 • S), lasting at least 2 days, were detected during this period, with vertical structures spanning two levels or more, i.e. multilevel events, during part of or for all their tracked life cycle. The number of cyclone events detected that corresponded at two levels is 347, at three levels 167, at four levels 99, at five levels or more 140. This does not mean that each event spanned all the levels by which it was categorized throughout its life cycle, rather that the category assigned represents the largest number of levels spanned at some stage in its evolution.
The validation is carried out in two steps. First, the overall seasonal tracking capability during the 2015 austral winter is compared to previous climatologies, in order to verify that the detected stacked cyclone events were geographically distributed, within reasonable limits of expected variability, in agreement with these climatologies. During 2015, a strong El Niño event took place with Multivariate El Niño/Southern Oscillation (ENSO) Index (MEI) values between 2 and 2.5 during the sampled period (https://www. esrl.noaa.gov/psd/enso/mei/table.html). The seasonal Southern Annular Mode (SAM) index reached a high positive value of 3.19 (http://www.nerc-bas.ac.uk/public/icd/gjma/ newsam.1957.2007.seas.txt), which may also imply changes to cyclone trajectories. The Indian Ocean Dipole index at this time also indicates a positive phase with comparatively large values (https://www.esrl.noaa.gov/psd/gcos_ wgsp/Timeseries/Data/dmi.long.data), while the Interdecadal Pacific Oscillation has weak negative values (https://www. esrl.noaa.gov/psd/data/timeseries/IPOTPI/). These variability indicators, known to impact SH weather systems (e.g. Pezza and Ambrizzi, 2003;Pezza et al., 2007;2008;Reboita et al., 2015), are considered here qualitatively since there may be nonlinear relationships with cyclone development and life cycle, as these studies point out. Their impacts on cyclone occurrence and development can only be evaluated through a fully developed climatology analysis.
The second validation step is focused on three different long-lived events over the SH during the winter (JJA) of 2015 spanning seven different pressure levels (925, 850, 700, 600, 500, 400 and 250 hPa). Two originate at subtropical/midlatitudes and one begins near Antarctica. These three events selected for SA validation were located at some point in their evolution in the vicinity of the Andes mountain range, which is a major orographic feature in the SH circulation (Inatsu and Hoskins, 2004). Following Gozzo et al. (2014) who included Geostationary Operational Environmental Satellite (GOES) satellite Infrared (IR) imagery in their case-study, both GOES IR, Water Vapour (WV) channels and ERA-Interim cloud-cover products are used here in order to validate the algorithm's capability and reliability.

Mean JJA 2015 state compared with climatology
The results for the 2015 austral winter are presented in Figures 4-8 and primarily but not exclusively refer to the SH cyclone climatology of Hoskins and Hodges (2005), which was based on TRACK and with the climatology of Lim and Simmonds (2007), which was based on the University of Melbourne pressure-based feature-tracking system (Murray and Simmonds, 1991a;1991b). Figure 4 (and Supporting Information Video S1) shows the wind speed at 300 hPa for JJA 2015 to highlight the mean position and variability of the jets during the period. The well-known SH jet structure can be seen with the spiralling of the subtropical jet into the polar jet, beginning over the eastern Pacific, over Chile and Argentina and into the South Atlantic, merging with the polar jet near the Greenwich meridian. Overall features agree with the Peña-Ortiz et al. (2013) jet-stream climatology using data from the NCEP/NCAR Reanalysis. The main differences are in the position of the transition of the subtropical to polar jet, which in 2015 occurs further south over southern South America and in a stronger polar jet close to Antarctica. These differences can be expected from the El Niño state and positive SAM during this particular austral winter. Figure 5 shows the cyclone event distribution during JJA 2015 for all multilevel cyclones detected. The largest cyclone density is located off the coast of Rio Negro and Chubut provinces, Argentina, over the South Atlantic, somewhat poleward of the axis of the poleward spiralling jet, in a region where it accelerates. Over the South Atlantic, close to the Greenwich meridian, two additional small high cyclone density areas can be seen, one equatorward and the other poleward of the jet axis. The other prominent cyclone density region is found close to Antarctica almost under the axis of the polar jet extending from the eastern South Atlantic across the southern Indian Ocean. Another small cyclone density region is located over the central Pacific under the subtropical jet, where it reaches its highest speed. Figure 6 specifically shows the cyclone density for "shallow" cyclones (spanning 2 and 3 levels) at low, mid-and upper troposphere levels, respectively, "intermediate" events  Figure 7) spanning lower or upper troposphere, respectively, and vertically "well-organized" cyclones ( Figure 8) spanning most if not all the levels considered in this study, building on Lim and Simmonds (2007) definitions and discussion. Even though only a single season is considered, which as noted is particularly impacted by a strong ENSO event, its mean state is comparable with the above-mentioned climatology studies, spanning at least 45 and 22 years respectively. Overall results found here agree with the broad features of cyclone occurrence over the SH winter, after taking into account the known consequences of the previously described state of the climate system.
Total shallow cyclone (events detected at 2 to 3 levels, Figure 6) densities show that, while present throughout the extratropical SH, the more frequent occurrence of such systems is on the Atlantic coast of southern South America. Densities also increase near the subtropical jet over the South Pacific and near the entrance of the subtropical jet over the South Atlantic. When considering the occurrence of shallow events in different height ranges, near-surface (925, 825 and 700 hPa surfaces, Figure 6b) densities maximize off the coast of Argentine Patagonia in the South Atlantic near Golfo Nuevo and Golfo San Matías between 40 • and 45 • S, just poleward of the spiralling jet. This is a prominent cyclogenetic region in the SH (Simmonds and Keay, 2000;Hoskins and Hodges, 2005). This region appears to extend into the central South Atlantic. Another region of prominent cyclone density would have been expected over the Rio de la Plata Basin in subtropical southeastern South America. However, present results agree with Pezza and Ambrizzi (2003) and Reboita et al. (2015) who noted (also using the University of Melbourne tracking code applied to MSLP alone) that surface cyclone occurrence is strengthened over the Patagonian Atlantic coast during El Niño events. The lack of a cyclogenetic region near the Rio de la Plata basin (e.g. Hoskins and Hodges, 2005) and such a strong feature along the coast of Patagonia could be linked to the widening of the Hadley Cell during El Niño events and the associated southward displacement of the subtropical jet, as seen in Figure 4. Another prominent region is found near Antarctica due south of Australia in the vicinity of Terre Adélie and George V Land, where the polar jet has weakened significantly ( Figure 4). Some shallow cyclone density secondary areas are found in the Indian Ocean extratropics, seemingly along the equatorward edge of the polar jet, which is strong in that longitude region. Note the central South Pacific region with reduced cyclone density, also in agreement with previous studies.
Mid-tropospheric (700, 600 and 500 hPa surfaces, Figure 6c) shallow systems primarily occur off the Atlantic coast of Patagonia, extending into the central South Atlantic, but with a narrower latitude extent than the surface shallow cyclone density. This region extends under the subtropical jet, where jet wind speed begins to increase prior to its spiralling transition to the polar jet, over the South Atlantic (Figure 4). Upper-level (500, 400 and 250 hPa, Figure 6d) shallow systems maximize to the west of South Africa between 20 • and 30 • S. These occur in the vicinity of the subtropical jet entry region west of Africa. A secondary maximum is also observed in the central South Pacific subtropics where maximum subtropical jet speeds are reached. As in the lower levels there is reduced activity in the central South Pacific. Upper-level shallow cyclone density (Figure 6d) off the coast of Patagonia is reduced with respect to the mid-and lower stratosphere. The first two density areas are in agreement with cut-off low (COL) occurrence in Pinheiro et al. (2017). However, as the authors note, the cyclone centres detected by TRACK in the upper troposphere may include pressure troughs as well as closed circulation systems at upper-tropospheric pressure levels, and hence further criteria need to be applied to confirm the occurrence of COLs. Once more, reduced shallow cyclone occurrence is observed in the central South Pacific at extratropical latitudes.
Intermediate cyclones are less frequent than shallow events (Figure 7a). They appear to occur more frequently over the Indian Ocean and over northern Chile, Argentina, Bolivia and Paraguay. While intermediate cyclones can be detected anywhere between subtropical and high latitudes, these systems (events detected at 4 to 5 levels) in the low to mid-troposphere (925, 850, 700, 600 and 500 hPa surfaces, Figure 7b) have larger densities only in the southern Indian Ocean between high and extratropical latitudes, primarily on the equatorward edge of the polar jet. Mid-to upper-level (700, 600, 500, 400 and 250 hPa surfaces, Figure 7c) intermediate events occur in the same areas. As was the case with shallow events, there is a reduced density region in the extratropical central South Pacific.
Well-organized cyclone events considered here are those that span most heights, i.e. 6 or 7 levels ( Figure 8) and are mostly found between 45 • and 60 • S, spanning an extended longitude range from the central South Atlantic into the southern Indian Ocean, as far as Australia. This result agrees with Lim and Simmonds (2007). These well-organized cyclone events have their largest density in the vicinity of the jet entry region and strongest section of the polar jet. A secondary area occurs where the polar jet weakens considerably. While such well-organized cyclone events may arise throughout the subtropical and extratropical SH, they appear to be essentially higher extratropical latitude features with very limited occurrence in the subtropics, at least during an El Niño year. Inspection of sea-ice satellite maps for the 2015 winter (not shown) appear to suggest that both the polar jet and well-organized cyclones are located near the edge of the sea-ice/open ocean transition.
Summarising, the detected cyclone systems cluster as expected in the vicinity of the subtropical and polar jets as well as near the spiralling transition of the jet over the South Atlantic. Most of the higher density clusters occur in regions where the upper-level jets accelerate or decelerate. Well-organized cyclones spanning the whole troposphere appear to occur also under the strongest winds in the polar jet. The SA results show that the cyclone systems spanning different height ranges tend to occur in different regions of the SH. The reasons for the varying spatial distribution of cyclones occurring at different height ranges, in particular for shallow and intermediate cyclones, is beyond the scope of the current study, requiring a further climatology study using a larger number of years.

Validation cases
Well-organized, long-lived events are considered here, given that they represent the full analysis of SA capability. The events analysed are merged bottom-to-top and top-to-bottom stacks. Even though the well-organized events span most or all of the seven pressure levels, SA results are displayed here at five levels only (925, 850, 700, 500 and 400 hPa) for the sake of clarity. In each case, a first set of images shows the stacking on a 3D display. The top image presents a map view of the stacking and the second image shows an oblique view of the same display, so that tracks at different heights can be viewed without overlapping. A second set of plots and images present both ERA-Interim cloud cover at 500 hPa and NOAA's GOES satellite imagery from the IR and WV imagery channels, respectively, as a four time-step sequence of evolution from the events. Each set also includes a close-up of the geopotential field at 300 hPa and 1,000-500 hPa thickness field, centred in the position of the RV anomaly at each time step. Finally, animation loops, for each of the three events, are provided as supplementary material to facilitate the visualization of cyclonic system evolution as captured by SA in ERA-Interim cloud cover and 300 hPa wind intensity together with GOES IR and WV images.

Event 1: 29 July-10 August 2015
The system for this validation case developed between 29 July and 10 August 2015, i.e. the RV anomaly life cycle lasted 12 days. It is based on 236 independent cyclone event detections by TRACK, combined by SA (Figures 9  and 10 This event initially appears to develop in the upper levels (400 and 500 hPa) and a few days later appears in the lower levels, moving east along the poleward edge of the subtropical jet (Figures 9 and 10a). At this point the RV anomaly stretches under the poleward edge of the subtropical jet between a confluent ridge to the north and a diffluent trough, aligned in a southeast-northwest axis, to the south. These features are visible both in the 1,000-500 hPa thickness and 300 hPa geopotential, though somewhat less intense in the latter. All these features suggest significant baroclinicity. Two days later this circulation structure continues as the system moves east. Note the existence of the semi-permanent anticyclone near tropical latitudes off the coast of Peru in this and subsequent time steps. By 7 August the trough sharpens as it becomes narrower, particularly in the 1,000-500 hPa thickness ( Figure 10c). Inspection of Figure 10c shows both in clouds (ERA-Interim and GOES), winds and RV the presence of a small cold front along the axis of the trough. Figure 10d shows a sharp trough aligned with the coast of Chile as the cyclone reaches maturity, as it moves close to the highest part of the Andes mountain range, with peaks above 6,000 m. The system rapidly veers south when it comes almost on top of the Andes, and it decays, following the jet. This southward motion is probably induced by a high-pressure system over central South America, visible in the WV image geopotential and thickness and animation. It follows the coast for about a couple of days until it vanishes. The RV anomaly showing a well-developed cyclone that is very organized through its vertical structure and the cyclone is clearly visible in the infrared and water vapour channel imagery. This event lacks a track at 250 hPa. It is assumed that the event did not develop at that altitude and evolved below that level. It may also be that either TRACK did not detect a track at that altitude or because the SA did not find during this period a track to link with the rest of the tracks in the stack. It is nevertheless tracked through six pressure levels, making it a distinct well-organized structure for a case-study.

Event 2: 27 June-11 July 2015
This validation case is tracked through all seven pressure levels and appears on all five levels displayed in Figures 11 and  12. This event is not always as distinct in GOES imagery as is the previous one, because of more dispersed cloud features. Hence it has a better signature in WV imagery than in the IR imagery, but it is still traceable and appears to respond differently to the presence of the Andes mountain range compared with the previous case. Initially the 300 hPa winds exhibit weak jet features (Figure 12a). It is first detected at 925 hPa on 27 June over the Pacific Ocean (32.95 • S, 225.17 • E), 6 h later at 850 hPa and the next day at all other levels (Figures 11 and 12). During the early stages it is associated with a weak diffluent trough followed by a somewhat stronger ridge in 300 hPa geopotential height and 1,000-500 hPa thickness fields (Figure 12a), but with no confluent ridge as in the previous event. It moves eastward and slightly equatorward, reaching southern South America in the vicinity of Puerto Montt (∼41 • S, Chile, Figure 12b) and the diffluent trough is somewhat stronger. During this stage the system evolves along the poleward edge of a weak polar jet (Supporting Information Video S3). At this somewhat higher latitude the average height of the Andes is limited to 2,500-3,000 m and thus the system is able to continue its development over and beyond the Andes. It crosses Argentina diagonally towards the Rio de la Plata (34 • S), reaching the South Atlantic near the border between Uruguay and Brazil. During this stage the system is tracked at 850 hPa and levels above. The GOES images and ERA-Interim cloud cover plots show a weak mature cyclone signature over the South Pacific, approximately 24 h before reaching the Andes. The cloud signature associated with the RV anomaly is weak after it crosses the continent, most probably because of low levels in atmospheric water content as it travels over the semi-arid Patagonia steppes particularly during winter. The RV anomaly couples over Buenos Aires with the poleward edge of the subtropical jet after becoming separated from the edge of a polar jet branch over Argentina. At this time the trough is sharper in the 1,000-500 hPa thickness map than at 300 hPa geopotential (Figure 12c). Also, a cloud signature begins to develop probably due to enhanced atmospheric water content as the anomaly travelled over the humid Pampas Region and the warm waters of the Rio de La Plata estuary. The cyclone cloud signature then reappears over the South Atlantic at a distance from the coast during the latter part of the event (Figure 12d). During this period the trough changes from a southeast-northwest slanted feature to a south to north feature suggesting cyclone maturity. The system evolves on the poleward side of the subtropical jet, under the influence of enhanced baroclinicity and probably atmospheric water vapour. The system vanishes over the mid-eastern Atlantic Ocean (20.98 • S, 35.73 • E). It is important to highlight the difference between the results of a stack analysis versus TRACK applied to a single level in the lower troposphere. In the latter case this particular event would have been detected as two independent events. In the first case, the STACKER algorithm is able to link them as an evolving single system which undergoes a first development over the South Pacific, almost has no signature in the drier atmosphere over Patagonia and then starts developing once more over the subtropical South Atlantic.

Event 3: 10 July-24 July 2015
This validation event is first detected close to the coast of eastern Antarctica. It is first tracked over the Antarctic July, and is present at the beginning at all levels with a few hours difference (Figures 13 and 14, Supporting Information Video S4). Of the three validation cases, it is the event which shows the most change in the levels of the tracks as the system evolves as well as the longest overall track since it virtually circumnavigates all the extratropical Southern Hemisphere. The lowermost 925 hPa track lasts a couple of days and vanishes by the time the system travels south of New Zealand (Figure 13b). By the time it reaches the mid-South Pacific, it is limited to 500 hPa and above, including 250 hPa (not shown). A probable cause for the initiation of this event could be found in the thermal difference between sea-ice and open waters near the initial tracking point, which is close to the edge of the sea-ice in the mean July 2015 satellite sea-ice retrievals. The system is tracked at and above 500 hPa as the system veers away from Antarctica, north of the Bellingshausen Sea, following the polar jet on its poleward side, as the latter moves north. As the RV anomaly nears southern South America (18 July, Figure 14a) a distinct cyclone cloud pattern can be seen both in IR and WV images. At this stage, the cyclone is surrounded by the polar jet. A diffluent trough can be observed both in the 1,000-500 hPa thickness and in 300 hPa geopotential height. Also note the ridge to the northwest of this trough near 35 • S. By 19 July ( Figure 14b) the system path reaches the continent and begins to cross the Andes north of 40 • S, where the Andes rapidly vary in height from 3,000 to 5,000 m towards the Equator. At this point the RV anomaly is no longer linked to the polar jet and reaches the subtropical jet. There appears to be a small front behind the system visible in the 300 hPa winds and a much stronger one preceding the system extending over Argentina and into the South Atlantic. Both have comparatively weak cloud structure, even more so the latter one. The cyclone system itself is only tenuously visible in the WV image. The tracks cross Argentina north of the Rio de la Plata near 32 • S and cross Uruguay and Brazil into the South Atlantic, veering to the southeast (21 July, Figure 14c). The cyclone structure begins to reappear as the system reaches the Atlantic Ocean, following the front. The system reaches 45-50 • S by 22 July (Figure 14d) and continues its life cycle at that latitude ( Figure 13). The subtropical jet resumes its more zonal behaviour near the Tropic of Capricorn. A high-pressure system develops to the northeast of the cyclone, with a confluent ridge. The cyclonic cloud pattern develops as the system moves over the South Atlantic, reaching maturity. The event finally decays 14 days after being first tracked, 24 July, very close to the original detection region (52.0 • S, 3.1 • W).

Additional comments about STACKER based on validation cases
A number of issues regarding SA come out from the above event analyses. Each validation case discussed above has a different, distinct evolution which SA was able to determine and show the differences. In the first place, note that the SA track points determined using RV are well co-located with respect to the evolving troughs in all three events. While this is linked to the specific tracking algorithm used, it shows that the tracks at various levels are well-correlated with relevant circulation features. Furthermore, the RV anomaly features found at the different levels, for a given event, are consistent with these observed regional/local circulation features. Figure 15 shows three such RV anomaly stacks at a given time, one for each of the three events discussed above. The figure shows that, be it a snapshot of a two-layer state of an originally multi-level event, or a multilevel event active at all levels, the RV anomalies in each case have a high degree of congruence with corresponding circulation features present at each level of an event as well as with each other at the different levels, highlighting how RV anomaly features vary with height within a cyclonic system. Furthermore, the comparison with both infrared and water vapour imagery which convey information of cloudiness and circulation on different height ranges validates the observed vertical structure found by SA for the cyclone events.
Inspection of the coordinates (latitude and longitude) at the various stacked levels for each event shows that as the initial RV anomalies evolve within the stack, the initial vertical slant July at 1800 UTC of the system (the RV anomalies at upper levels lag behind those at lower levels) then diminishes as cyclonic features develop and mature. The differences in latitude and longitude between upper-and lower-level tracks become very small, in agreement with the baroclinic instability theory for cyclogenesis. Furthermore, the inspection of the troughs observed in 1,000-500 hPa thickness and 300 hPa geopotential shows that as the RV anomalies evolve into distinct cyclonic features, their inclinations with respect to meridians decrease until the troughs finally becomes parallel to them, when the systems reach full development.
Consequently, the SA is correctly combining tracks at different levels corresponding to the same cyclone event and provides accurate 3D views of the evolving systems.

CONCLUSION
A 4D-STACKER algorithm has been developed, following dynamic programming criteria, to track extratropical RV anomaly 3D structures over time with an objective systematic approach. The algorithm is based on the TRACK algorithm which provides the raw tracks at single levels, which must then be combined objectively to determine the 3D structure and its evolution. Because of the extremely large number of possible single-level track combinations, an optimal solution is reached according to a defined criterion necessary to reduce options and better define the track "stacks" representing a given cyclone event. The SA is progressive, beginning with two initial levels and then adding additional levels to the stack. This allows an iterative stacking approach, adding one level at a time, which furthermore allows for a bottom-up or top-down approach. Furthermore, bottom-up and top-down results can be objectively merged to obtain a better description of the detected events. This leads to an improved 4D view of the RV anomaly events, and hence of the structure of cyclones necessary for the advancement of atmospheric system dynamics and dynamic meteorology.
SA is validated through its application to a winter season (JJA 2015) over the extratropical Southern Hemisphere, and in particular southern South America, where the Andes mountain range provides an orographic barrier that strongly impacts on the SH atmospheric circulation (Inatsu and Hoskins, 2004). SA validation is specifically carried out using cyclone tracks and features. Their behaviour is impacted by the Andes in different ways depending on latitude, giving valuable information on the capability of SA to capture such features.
After taking into account specific issues pertaining to currently known impacts of major variability processes such as ENSO on the SH extratropical cyclone development, SA results are in very good agreement with current cyclone climatologies mostly based on single-level studies and limited to the lower troposphere. The behaviour of the cyclones and their tracks during the sampled period verifies the known dynamic behaviour of cyclones and their link with upper-atmospheric jets. Cyclone densities for all events agree well with the climatologies, and the specific analyses for the shallow, intermediate and well-developed cyclones show that the different vertical organizations each have preferred zones of occurrence.
The three case-studies, based on long-term well-developed systems, highlight the impact of the Andes depending on where the system reaches the mountain range. The Altas Cumbres region of the Andes (6,000 m and above) appears to disrupt the system. Keable et al. (2002) and Hoskins and Hodges (2005) have identified this as a 500 hPa cyclolysis region. The lower Patagonian Andes have a reduced impact limiting the development of the well-developed cyclone structure but once the systems cross the mountain range and enter the South Atlantic the development of the cyclone structures starts once more. Such development appears to take place near well-known cyclogenetic regions. These validation studies highlight the importance of looking at the 3D temporal evolution of individual cyclones and other vorticity anomaly events in order to adequately understand their behaviour.
Consequently, the STACKER Algorithm, using TRACK as the feature-tracking algorithm, demonstrates the applicability of dynamic programming for both 4D dynamic meteorology and dynamic climatology studies, as evidenced here by the two-stage validation. Understanding how complex 3D structures such as cyclones evolve in dynamically, spatially varying flows, and the classification of their main multilevel aspects and behaviours, remains a significant challenge. Furthermore, understanding how climate variability and climate change impact cyclone and severe weather events and their distribution also requires a 4D perspective. SA is therefore a valuable tool for studies of vortex dynamics. A 4D-STACKER-based SH dynamic climatology is currently under way. relative vorticity (or other anomalies in T, Z or P) by TRACK during the preliminary stage of the analysis. Detections are assumed to be meteorological events, e.g. cyclones, anticyclones or other significant perturbations in the atmosphere, at a given instant or time frame. Detections in TRACK are not directly related to each other until the tracking optimization process is completed. Points are associated into Tracks by the tracking algorithm, e.g. TRACK.

A2 TRACKS
Tracks are the main object obtained by the tracking algorithm. They consist of a series of points, from consecutive time steps associated according to the criteria that weather features cannot change their position, direction, velocity or intensity abruptly. It is safe to assume that a given Track represents at least a fraction of a featurés translational behaviour during its life cycle. A set of points defining a Track will always be consecutive.

A3 STACKS
Stacks are the main output of the STACKER algorithm. Stacks consist of a series of associated Tracks on several levels. These are four-dimensional objects, in a combination of discrete values within the longitude, latitude, height and time ranges chosen for the study (cf. Figure 3). Stacks are always composed of Tracks stacked from consecutive levels (levels may correspond to pressure, model, PV or theta levels), i.e. the stacked detected event must be continuous in the vertical axis. It is possible to find several non-overlapping Tracks associated on the same level, because of the Point-to-Point linking approach, e.g. Stack 4 at the 850 hPa level in Figure 3. Tracks on a given Stack are directly related only to Tracks on the first-lower and first-upper level and are indirectly related to the rest of the Tracks in the Stack.

A4 SIMILITUDE FUNCTION AND THRESHOLD VALUE
The most important feature of the method is the design of the objective cost function and the selection of threshold values. The criterion to determine the likelihood of two Tracks being associated is defined in the Similitude Function (SF). It assigns a value to each optimization variable, or Link, and penalizes, as a cost, dissimilar behaviour between detections, e.g. distance between Points, significant deviations in velocity or directions between Tracks, differences in size or intensity of the events, etc. This "Similitude Cost" must have a Threshold which states that any Link between two Points that have a high SF value is not relevant to the problem. If the two Points in a Link belong to two Tracks such that, at some point, two of its Points are possible links according to the SF criteria, then the Link should be considered relevant for the study (Figure 3). The threshold will act both as a filter to avoid the build-up of a massive numerical problem to solve, and provide a "default" SF value for the solver to choose if there are no better options to associate a Point with another one.

A5 LINKS
Links represent the hypothesis that two single and simultaneous detections from two consecutive pressure levels correspond to the same meteorological event. These are the association between any lower Point and any simultaneous upper Point. Each link is assigned a cost value according to the Similitude Cost Function (SCF), depending on the classification of the Link. There are as many optimization variables as there are Links created in the analysis. The solving process will decide which Links are more likely to couple two detections of the same meteorological event at adjacent levels. Since Links are rather complex compared to the other objects in the algorithm, they are classified according to their behaviour as True, False and Ghost Links.

A6 SIMILITUDE COST FUNCTION
The Similitude Cost Function (SCF) sets the priority value on Links to be activated later in the solving process and is derived from the SF. Its value is assigned to a Link and later used as a coefficient in the Objective Function. Since the objective function must be minimized, lower cost values will be prioritized over others. Hence the importance of its criterion and Threshold value. Its desired behaviour is to prioritise True Links over the rest, penalize False Links, which exceed the present Threshold and present Ghost Links as "default" option to the solver. More specifically:

A.6.1 Threshold value
In this study, Points are limited to longitude and latitude, in order to calculate distance. For more complex studies the core set of parameters can include other variables, e.g. altitude, wind speed, pressure, temperature, humidity, etc. The Threshold value is arbitrary and is specifically chosen for each study according to its aims, precision required and processing resources available. Here, following Lim and Simmonds (2007), the SF is equal to the "distance" between a lower and an upper Point, so the Threshold value can be understood as the radius of a circular surface centred at the projection of the Lower Point's position on the upper level. "Distance" is defined as the arc length between the projections of the Points on a sphere with Earth́s mean radius.

A.6.2 True links
These are the essence of what Links represent. A True Link is generated between two simultaneous Points when its SF value does not exceed the Threshold value. It means that the relation between these two Points is relevant for the algorithm according to the current study criteria. Its Cost value is assigned directly from the Similitude Function and it will be smaller than the Threshold value.

A.6.3 False links
False Links complement True Links. False Links are Links between two simultaneous Points where their SF value exceeds the Threshold value. The algorithm generates them because each Point belongs to a Stack that, at some time step, is related by a True Link with the other Point's Stack.
Since the Stacker Algorithm is not allowed to associate single simultaneous Points without all of the other simultaneous Points between two Stacks, the algorithm False links complement what could not be True Linked. Its Cost value is equally assigned by the SF and it will be larger than the Threshold value, thus penalizing associations with a large number of these Links and prioritizing associations with larger numbers of True Links.

A.6.4 Ghost links
Ghost Links relate two simultaneous Points, but one of them does not represent a local maximum/minimum detection from the Tracking Algorithm. These objects are artificially generated in order to give the linear programming solver a default alternative to choose when there are no better options to associate a Point with. Ghost Points are essentially a projection of a lower Point in the upper level in a given stage of the algorithm and do not belong to any Stack. The Stacker Algorithm specifically assigns the Threshold value to its Cost value in this case.

A.6.5 Binary linear programming
The STACKER Algorithm is a progressive optimization method. The whole problem is modelled so that it can be solved through Dynamic Programming by disassembling it into smaller problems. To achieve its final solution, it iterates through a series of local stages which only includes the association of Points from two consecutive pressure levels. Each of these local problems are proposed as Linear Programming problems to guarantee the existence of an optimal solution. These problems can be solved by a specific case of Linear Integer Programming, called Binary Programming, in which the optimization variables can only take two possible values, either 1 or 0, i.e. activated or deactivated. This way the solver is proposed with several options to activate or deactivate (Links) according to the criteria of minimizing a Similitude Cost Function, respecting the restrictions that constraints provide.

A.6.6 Constraints
Constraints are, by far, the most complex and delicate features of the algorithm because of their progressively increasing size. These constraints condition the solver's behaviour in order to achieve the desired results, each addressing a specific need:

A.6.7 Must Be Together Constraint (MBT)
This constraint is responsible for assuring that Stacks will not be disassembled into smaller groups of Points by the solver.
Since the aim of the algorithm is a Stack-to-Stack association using Point-to-Point links, if a single Link between an upper and lower Stack has been activated by the solver, all Links between those two same Stacks must be activated too.

A.6.8 Must Not Be Together Constraint (MNBT)
This set of constraints prevents the solver from associating Stacks that overlap with each other at some point (Figure 2), i.e. two different Stacks that contain at least a single Point on the same surface and time step, must be different events, and thus, must not be combined into the same Stack.

A.6.9 Do Not Repeat Lower Tracks (DRLT)
This set of constraints forces the solver to activate at least one Link for every lower Point and makes it impossible for the solver to associate the same lower Point to more than one upper Point. In this case, for every set of Links containing the same lower Point, only one must be activated at the same time. The activation of a True or False Link represents the association between two Stacks. The activation of a Ghost Link indicates that the solver took the "default option" and any other Links available for its lower Point were not optimal for the solution of the local problem.

A.6.10 Do Not Repeat Higher Tracks (DRHT)
Like the previous set of constraints, DRHT makes it impossible for the solver to associate the same upper Point to more than one lower Point. The only difference lies in the fact that it will not force the solver to activate at least one Link for every upper Point.

A.6.11 Binary solution parsing
Once the Linear Programming process is complete, the solver returns all activated and deactivated Links. Deactivated Links are irrelevant to the parsing, as well as active Ghost Links and have no practical meaning. Active True and False Links represent the most coherent associations between lower and upper Stacks according to the SF criteria. Parsing must take all of the associated Stacks' Points and build a new set of Stacks on as many levels as the previous two sets combined. Point-to-Point linkage allows the association of any number of Stacks together: • A Stack may have no coherent associations, and hence no changes are required. • Any number of lower Stacks can be associated with any number of higher Stacks, building a new Stack with as many levels as the previous Stacks combined had and all of their Points.
New Stacks will become lower (upper) Stacks on the next iteration of the problem or will become the final solution if there is no further stacking possible. If the process is repeated in the opposite vertical direction (Top-to-Bottom or Bottom-to-Top) the Stacker Algorithm will provide a similar but somewhat different solution. Merging Top-to-Bottom and Bottom-to-Top stacks combines all detected features into a single stack for subsequent analysis.

A.6.12 Stack concatenation
In this study, the Linear Programming package from the R programming language was used to implement the algorithm. Due to computational limitations, the size of the problem becomes unsolvable for large amounts of variables (over 10,000) and a Stack concatenator process is proposed. Concatenation consists of identifying the same meteorological event through several multi-level Stacks obtained in short samples and combining all of the sections into a single, longer Stack. To this end, we used an overlap of 5 days between datasets and assumed that Stacks that shared at least one detection (Point) during the overlapping period must be the same event and thus must be combined. This is possible since in principle the analysis of fields by tracking algorithms does not suffer from time sampling edge effects. The concatenation process is an approximation and may produce some overlapping and occasional discrepancies in the combined Stacks. Since this impacts less than 0.005% of the whole Stack set, it does not compromise the coherence of the algorithm's results.