A simple kinematic model for the Lagrangian description of relevant nonlinear processes in the stratospheric polar vortex. Nonlinear Processes (2), 265-278.

. In this work, we study the Lagrangian footprint of the planetary waves present in the Southern Hemisphere stratosphere during the exceptional sudden Stratospheric warming event that took place during September 2002. Our focus is on constructing a simple kinematic model that re-tains the fundamental mechanisms responsible for complex ﬂuid parcel evolution, during the polar vortex breakdown and its previous stages. The construction of the kinematic model is guided by the Fourier decomposition of the geopotential ﬁeld. The study of Lagrangian transport phenomena in the ERA-Interim reanalysis data highlights hyperbolic trajectories, and these trajectories are Lagrangian objects that are the kinematic mechanism for the observed ﬁlamentation phenomena. Our analysis shows that the breaking and splitting of the polar vortex is justiﬁed in our model by the sudden growth of a planetary wave and the decay of the axisymmetric


Introduction
The availability of high-resolution and high-quality reanalysis data sets provides us with a powerful tool for obtaining a detailed view of the space-time evolution of the stratospheric polar night vortex (SPV), which has implications for the geophysical fluid dynamics of the entire Earth. The complexity of such a detailed view, however, makes it difficult to extract the physical mechanisms underlying notable transport features in the observed behaviour. The goal of this work is to gain new insights into the fundamental mechanisms responsible for complex fluid parcel evolution, since these lie at the heart of our understanding of the dynamics and chemistry of the stratosphere. To this end, we extract, directly from the data, a simple model with strippeddown dynamics in order to directly probe, in a controlled and systematic manner, the physical mechanisms responsible for the key observed transport features of the SPV. Models of this kind, termed "kinematic models", have provided a simple approach for studying Lagrangian transport and exchange associated with flow structures such as meandering jets and travelling waves (Bower, 1991;Samelson, 1992;Malhotra and Wiggins, 1998;Samelson and Wiggins, 2006). Other works have used analytical kinematic models to illustrate phenomena in planetary atmospheres (e.g. Rypina et al., 2007;Morales-Juberías et al., 2015). In the present paper, we focus on SPV transport processes associated with filamentation and vortex breaking, of which the dynamical structure is not fully understood.
The importance of an increased understanding of the SPV was dramatically demonstrated by the intense research effort that followed the discovery of the "Antarctic ozone hole" phenomenon in the 1970s (Chubachi, 1984;Farman et al., 1985;Solomon, 1988). In the following decades, during which monitoring of ozone in atmospheric columns above Antarctica showed little interannual variability, in situ measurements corroborated by satellite data revealed that ozone was systematically decreasing in the Antarctic lower stratosphere during the southern spring season. Whilst this was immediately associated with the simultaneous increase in atmospheric pollution by anthropogenic activities, several key questions arose (Solomon, 1999): 1. Why does this occur over Antarctica and not over the Arctic (since pollution sources are stronger in the Northern than in the Southern Hemisphere)?
2. Why does this occur in the spring season?
The research demonstrated that, indeed, increased atmospheric pollution was to be blamed for the ozone depletion and identified the participating substances and special mechanisms. The research also demonstrated that the unique atmospheric conditions above Antarctica were responsible for the geographic preference for ozone destruction. In particular, it was shown that the strong circumpolar and westerly SPV characteristic of the southern winter and spring stratosphere contributes to the isolation of the cold polar region, setting up a favourable environment for the special chemistry to act. The new knowledge led to the formulation of international agreements that resulted in a negative answer to question (3) above. The analysis of transport of fluid parcels outside the region isolated by the SPV also showed strong stirring and mixing of the flow. In this "surf zone" (McIntyre and Palmer, 1984), air parcels can travel long distances away from the SPV in an environment where contours of long-lived tracers, such as potential vorticity, can stretch and form complex patterns. In this region, Rossby wave breaking (RWB) is associated with irreversible deformation that pulls material filaments of the outer edge of the SPV and enhances mixing with the exterior flow (McIntyre and Palmer, 1983, 1984. Such a process makes the SPV edge a barrier to horizontal transport of air parcels (Juckes and McIntyre, 1987) while continuously eroding and regenerating the SPV edge by filamentation (Bowman, 1993). Polvani and Plumb (1992) and Nakamura and Plumb (1994) examined, in an idealized setting, the way in which Rossby waves break and eject SPV material outward. The latter conceived a similar setting in which Rossby waves also break inward. Dynamical systems theory provides valuable insights into the transport processes described in the previous paragraph. Tools of the theory include the geometrical structures, referred to as hyperbolic trajectories (HTs), their stable and unstable manifolds, and their intersection in homoclinic and heteroclinic trajectories that provide the theoretical and computational basis for describing the filamentation process. A challenge in the application of these concepts to realistic geophysical flows is that while the structures mentioned are defined for infinite-time autonomous or periodic systems, geophysical flows are typically defined as finite-time data sets and are not periodic. Mancho et al. (2006b) addressed this challenge for realistic ocean flows by identifying special hyperbolic trajectories in the finite data set, called distinguished hyperbolic trajectories (DHTs), and by computing stable and unstable manifolds as curves advected by the velocity field. A pioneering effort for identifying HTs for the stratosphere was due to Bowman (1993). McIntyre and Palmer (1983), Bowman (1996) and de la Cámara et al. (2013) suggested that HTs are responsible for the cat-eye structures associated with planetary wave breaking at the critical levels, i.e. where the wave phase speed matches the background velocity (Stewartson, 1977;Warn and Warn, 1978). HTs are at the locations where the cat-eye structures meet. Perturbation of the cat eyes results in irreversible deformation of material contours, signifying Rossby wave breaking. de la Cámara et al. (2013) and Guha et al. (2016) identified HTs both within and outside the SPV, thus suggesting that Rossby wave breaking can occur in either of those regions. The former authors worked with reanalysis data, while Guha et al. (2016) used a dynamical model based on the shallow-water equations in which the perturbing waves are produced in a controlled manner. Therefore, HTs are essential features for tracer mixing both outside and inside the vortex, and for occasional air crossings of the vortex edge.
We focus on the SPV behaviour during the major stratospheric sudden warming that occurred in the southern stratosphere during September 2002. In this unusual event, the SPV broke down in the middle stratosphere (Mechoso et al., 1988;Varotsos, 2002Varotsos, , 2003Varotsos, , 2004Allen et al., 2003;Konopka et al., 2005;Esler and Scott, 2005;Manney et al., 2006;Charlton et al., 2006;Taguchi, 2014). We begin by identifying key Lagrangian features of the flow in reanalysis data fields. Next, we build a kinematic model of the event that emulates the behaviour of planetary waves observed in the data. We show that our model produces strikingly similar transport features to those found in the reanalysis data, confirming the key role played by the HTs during vortex filamentation and breakdown.
The structure of the paper is as follows. Section 2 describes the data and methods we used. Section 3 describes the planetary waves in the reanalysis data in the year 2002 in the stratosphere at selected pressure levels (10 hPa). We relate these to filamentation phenomena and the polar vortex breakdown that occurred in that year. Section 4 reproduces the findings obtained with our analytical kinematical model, confirming the role played by the HTs in the 2002 vortex filamentation and breakdown. Section 5 discusses the consistency of the kinematic model as representative of atmospheric flows that conserve potential vorticity. Finally, in Sect. 6, we present the conclusions.

ERA-Interim reanalysis data
To achieve a realistic representation of the atmospheric transport processes, it is crucial to use a reliable and high-quality data set. We use in this work the ERA-Interim reanalysis data set produced by a weather forecast assimilation system developed by the European Centre for Medium-Range Weather Forecasts (ECMWF; Simmons et al., 2007). de la Cámara et al. (2013) obtained encouraging results on the suitability of the ERA-Interim data set for Lagrangian studies of stratospheric motions in their comparison of parcel trajectories on the 475 K isentropic surface (around 20 km) using this data set and the trajectories of super-pressure balloons released from Antarctica by the VORCORE project during the spring of 2005 (Rabier et al., 2010).
The ERA-Interim data set that we selected for this study is available four times daily (00:00, 06:00, 12:00 and 18:00 UTC), with a horizontal resolution of 1 • × 1 • in longitude and latitude and 60σ levels in the vertical from 1000 to 0.1 hPa. The data cover the period from 1979 to the present day (Dee et al., 2011) and they can be downloaded from http://apps.ecmwf.int/datasets/data/ interim-full-daily/levtype=sfc/. In particular, we will use the data for the geopotential height on surfaces of constant pressure and wind fields on isentropic surfaces for the period August-September 2002.
The geopotential height Z on constant pressure surfaces p is defined as the normalization to g 0 = 9.80665 ms −2 (standard gravity at mean sea level) of the gravitational potential energy per unit mass at an elevation s (over the Earth's surface) and has the form where g is the acceleration due to gravity, λ is longitude, φ is latitude and z is the geometric height (Holton, 2004). In the quasi-geostrophic approximation, the geopotential height is proportional to the streamfunction of the geostrophic flow (Holton, 2004). For the analysis of planetary waves, we apply a zonal Fourier decomposition to the geopotential height field on the 10 hPa pressure level (approximately 850 K potential temperature). The zonal wave decomposition yields (2) The mean flow is defined as and the different modes Z k with wave number k ≥ 1 have the sinusoidal description: where λ ∈ [0, 2π ) is longitude, φ ∈ [−π/2, π/2] is latitude, B k is the amplitude of the wave and ϕ k its phase. During the warming event that occurred in the southern stratosphere during September 2002, the flow was dominated by the contributions of the mean flow and the two longest planetary waves (Z 1 and Z 2 ; Krüger et al., 2005)

Lagrangian descriptors
Dynamical systems theory provides a qualitative description of the evolution of particle trajectories by means of geometrical objects that partition the phase space (the atmosphere in our case) into regions in which the system shows distinct dynamical behaviours. These geometrical structures act as material barriers to fluid parcels and are closely related to flow regions known as hyperbolic, where rapid contraction and expansion takes place. Several Lagrangian techniques have been developed in order to detect such structures in geophysical fluids. This is challenging because, while classical dynamical systems theory is defined for infinite-time autonomous or periodic systems, in geophysical contexts, the velocity fields are generally time dependent, aperiodic in time and defined over a finite discrete space-time domain. Among others, techniques developed are finite-size Lyapunov exponents (FSLEs) (Aurell et al., 1997) and finitetime Lyapunov exponents (FTLEs) (Haller, 2000;Haller and Yuan, 2000;Haller, 2001;Shadden et al., 2005). Other techniques include DHTs (Ide et al., 2002;Ju et al., 2003) and the direct calculation of manifolds as material surfaces (Mancho et al., 2003(Mancho et al., , 2004(Mancho et al., , 2006b, the geodesic theory of Lagrangian coherent structures (LCS) (Haller and Beron-Vera, 2012) and the variational theory of LCS (Farazmand and Haller, 2012), etc. Our choice in this work will be the use of the Lagrangian descriptor (LD) function M introduced by Madrid and Mancho (2009) and Mendoza and Mancho (2010). The function M has been applied in a variety of geophysical contexts. For example, in the ocean, it has been used to analyse the structure of the Kuroshio current (Mendoza and Mancho, 2012), to discuss the performance of different oceanic data sets (Mendoza et al., 2014), to analyse and develop search and rescue strategies at sea (García-Garrido et al., 2015) and to efficiently manage in real time the environmental impact of marine oil spills (García-Garrido et al., 2016). In the field of atmospheric sciences, M has been used to study transport processes across the southern SPV and RWB de la Cámara et al., 2012Cámara et al., , 2013Smith and McDonald, 2014;Guha et al., 2016 and to investigate the Northern Hemisphere major stratospheric final warming in 2016 (Manney and Lawrence, 2016). The dynamical system that governs the atmospheric flow is given bẏ where x (t; x 0 ) represents the trajectory of a parcel that, at time t 0 , is at position x 0 , and v is the wind velocity field. Since our interest is in the timescale of stratospheric sudden warmings (∼ 10 days), we can assume to a good approximation that the fluid parcels evolve adiabatically. Therefore, trajectories are constrained to surfaces of constant specific potential temperature (isentropic surfaces). We will concentrate on the 850 K surface, which is in the middle stratosphere and approximately corresponds to the 10 hPa level. In Sect. 3, we expand on the reasons for this choice.
To compute fluid parcels trajectories, it is necessary to integrate Eq. (5). As the velocity field is provided on a discrete spatiotemporal grid, the first issue to deal with is that of interpolation. We apply bicubic interpolation in space and thirdorder Lagrangian polynomials in time (see Mancho et al., 2006a for details). Moreover, for the time evolution, we have used an adaptive Cash-Karp Runge-Kutta method. It is important to remark that, as done in de la Cámara et al. (2012) for the computation of particle trajectories, we use Cartesian coordinates in order to avoid the singularity problem arising at the poles from the description of the Earth's system in spherical coordinates. For our Lagrangian diagnostic, we use the M function defined as follows: where · stands for the modulus of the velocity vector. At a given time t 0 , the function M(x 0 , t 0 , τ ) measures the arc length traced by the trajectory starting at x 0 = x(t 0 ) as it evolves forward and backward in time for a time interval τ . Sharp changes of M values (what we call singular features of M) occur for sufficiently large τ , for very close initial conditions and highlight stable and unstable manifolds. Mancho (2010, 2012) have performed systematic numerical computations of invariant manifolds and found that they are aligned with singular features of M. They also provide examples in geophysical flows where manifolds are defined in a constructive way. Invariant manifolds are mathematical objects classically defined for infinite time intervals. The unstable (stable) manifold of a hyperbolic fixed point or periodic trajectory is formed by the set of trajectories that in negative (positive) infinity time approach these special trajectories. In geophysical contexts, this definition is not realizable, because only finite-time aperiodic data sets are possible. Nevertheless, manifolds can still be defined constructively with the following procedure. At the beginning time, these curves are approximated by segments with short length, aligned with the stable and unstable subspaces of the DHT identified with algorithms described in Ide et al. (2002) and Madrid and Mancho (2009). This starting step aims to build a finite-time version of the asymptotic property of manifolds. Next, segments are advected forward and backward in time by the velocity field. Due to the large expansion and contraction rates in the neighbourhood of the DHT, the curves grow rapidly in forward and backward time and specific issues are addressed by the procedure described in Mancho et al. (2003) and Mancho et al. (2004). The procedure provides curves, manifolds, that by construction are barriers to transport in geophysical flows. In this way, since manifolds are aligned with singular features of M, the latter belong to invariant curves of system Eq. (5), and therefore their crossing points are indeed trajectories of system Eq. (5). The capability of LDs, in general, and M in particular, to reveal invariant manifolds was analysed in detail in Mancho et al. (2013). Lopesino et al. (2015) and Lopesino et al. (2017) have discussed, in discrete and continuous-time dynamical systems, respectively, a theoretical framework for some particular versions of LDs in specific examples.
The consistency between the output field of Eq. (6) and FTLE ridges has been discussed in some references (see Mendoza and Mancho, 2010;de la Cámara et al., 2012;Mancho et al., 2013). The integral expression in Eq. (6) can be split in two terms: one for forward time and the other for backward time integration. Explicit calculations discussed in Mancho et al. (2013) for a linear saddle, show that singular features of the first term are aligned with the stable manifolds while those for the backward time integration are aligned with the unstable manifolds. This is similar to what is obtained with FTLEs that highlight stable and unstable manifolds, respectively, for forward and backward time integration intervals. The fact that we choose to add both fields is advantageous for highlighting hyperbolic trajectories at the crossing points of the singular features.
As an example relevant to the case that motivates the present study, we show in Fig. 1 the evaluation of M over the Southern Hemisphere using τ = 15 on the 850 K isentropic level for 5 August 2002. The representation shows a stereographic projection (see Snyder, 1987) in which the SPV is clearly visible by the bright yellow colour as well as the filamentation phenomena ejecting material both from the outer and inner parts of the jet. These filaments are related to the presence of hyperbolic trajectories highlighted in the figure. The fact that these saddle points of the LD field are hyperbolic trajectories of system Eq. (5) is numerically supported. To this end, de la Cámara et al. (2013) show that, for similar ERA-Interim fields, these points belong to the intersection of stable and unstable manifolds highlighted by the singular features of the field (see their Fig. 2). In what follows, all figures showing M were computed with τ = 15. This choice of τ is made based on the fact that diabatic heating/cooling processes in the extratropical stratosphere generally have longer timescales than those of horizontal advection. Hence, air parcels move on two-dimensional isentropic surfaces to a good approximation (they stay within 850 K for 30 days) (Plumb, 2007). Moreover, diabatic heating rates in the Antarctic mid-stratosphere are on the order of 0.5 K day −1 , although uncertainties in this magnitude remain large (Fueglistaler et al., 2009). During the time interval of our calculations of isentropic trajectories (τ = 15 days, i.e. time period of 30 days), the material surface would experience an increase of potential temperature of around 15 K. Nevertheless, calculations of M using wind fields at 850 K and 700 K (not shown) produce qualitatively similar results. This suggests that horizontal motions of the parcels will be affected by similar geometric structures at those isentropic levels and that the isentropic approach is justified in our problem.

Data analysis
As we indicated in the previous section, in order to characterize the planetary waves that propagate in the stratosphere, we carry out a Fourier decomposition of the geopotential height. In Fig. 2, we show the axisymmetric mean flow together with waves 1 and 2 in the geopotential field for 22 September 2002 on the 10 hPa pressure surface. The time evolution of these waves is described in the movies S1-S4 in the Supplement. Movies S1-S3 show components 0, 1 and 2 separately for the time period of interest, while S4 shows the superposition of these three waves. It is important to reiterate that, since the geopotential provides a good approximation of the streamfunction of the large-scale flow in the extratropical regions, its analysis will provide us with guidance on the building of the simple kinematic model presented in the next section. On the 10 hPa pressure level, the winter SPV in the Southern Hemisphere can be broadly defined as a circumpolar westerly jet. Figure 3a and b illustrate the evolution of the circulation during August-September 2002. We can clearly see the gradual deceleration of the SPV and the abrupt change in direction from westerly to easterly velocities at high latitudes that occurred on 22 September. This was a unique major sudden stratospheric warming (SSW) in the southern stratosphere. Planetary waves in the southern stratosphere were very active during the period where the 2002 SSW developed. Figure 3c presents a time series of the ratio between the amplitudes of waves 1 and 2. Increased wave 1 amplitude results in a displacement of the SPV vortex from a circumpolar configuration, while increased wave 2 amplitude results in a stretching the SPV in one direction and contraction (or "pinching") in the orthogonal direction. According to Fig. 3c, the amplitude of wave 1 was generally larger than that of wave 2 during the entire period, confirming the major role of this wave. Finally, Fig. 3d displays the variations in time of the ridges of wave 1 and wave 2. Note that wave 1 is quasi-stationary, while wave 2 propagates eastward, as is typical in the southern stratosphere during early spring (Manney et al., 1991;Quintanar and Mechoso, 1995).
The contribution of these different waves to the evolution of the SPV and their transport implications are clearly observed in movie S5. A regime giving rise to the stretching of material lines and the appearance of hyperbolic regions and the associated filamentation processes is observed. These filamentous structures and HTs are clearly highlighted by the application of LDs to the wind fields, as shown in Figs. 1 and 4. Filamentation phenomena occur both inside and outside the vortex, where the outer filamentous structures play the role of eroding the jet material barrier. Also, the presence of HTs in the flow (see captions of Figs. 1 and 4) indicates regions subjected to intense deformation and mixing (see Ottino, 1989). We emphasize that HTs appear both inside and outside the SPV. Finally, the breakup of the SPV on 24 September 2002 depicted in Fig. 4b (see also movie S5) occurs when manifolds associated with an HT that forms within the SPV connect the interior and the exterior of the jet, allowing for the interchange of parcels through the barrier. The pinching of the SPV takes place off the pole because Z 1 has large amplitudes in the days preceding the breakup. As we approach 24 September, Z 2 becomes the same order as Z 0 , and the jet elongates and flattens. At this point, the mean flow reversal is crucial for completing the pinching process

The kinematic model
Kinematic models have a long history in the geophysical fluid dynamics community. They allow for a detailed parametric study of the influence of identified flow structures on transport and exchange of fluid parcels. All early studies utilizing the dynamical systems approach for understanding Lagrangian transport and exchange associated with flow structures, such as meandering jets and travelling waves, have employed kinematic models (see Samelson and Wiggins, 2006).
Continuing in this spirit, we propose a kinematic model that allows us to identify, in a controlled fashion, the characteristics of the distinct propagating waves that are responsible for the different Lagrangian features observed in the SPV. Our kinematic model is inspired by the Fourier component decomposition of the geopotential extracted from the ERA-Interim data as discussed in the previous section. The analysis of data from August and September 2002 shows a mean axisymmetric flow, disturbed mainly by waves with planetary wave numbers 1 and 2 whose amplitudes and phase speeds vary in a time-dependent fashion. Therefore, we propose a kinematic model in the form of a streamfunction given by where ε 0 , ε 1 , ε 1 are the perturbation parameters, which we will refer to as amplitudes, and i are the Fourier components, along the azimuthal direction with wave numbers i = 0, 1, 2 respectively, which we describe next. We will work in a plane (x, y) that is the orthographic projection of the Southern Hemisphere onto the equatorial plane (Snyder, 1987). For simplicity, and in order to highlight the periodicity along the azimuthal direction, the components of the streamfunction are given in terms of polar coordinates satisfying x = r cos(λ) and y = r sin(λ), where the azimuthal direction λ is related to the geographical longitude and r is related to the geographical latitude.
The particular forms of 0 , 1 and 2 are inspired by the Fourier decomposition of the geopotential field shown in Fig. 2 for the 10 hPa pressure level on 22 September 2002. Starting with the mean zonal velocity, we will assume a jet with the following expression: Therefore, v λ = 0 only at r = 0 and r = a, and the velocity decreases exponentially away from the pole. Changing the values of a will allow us to consider variations in the position The other streamfunction components are 1 = −r 2 e −r 2 sin(λ) (10) and 2 = (r/d) 2 e −r 2 /d sin(2λ + ω 2 t + π/4), where d and w 2 are also tunable constants, and the phase π/4 was added so that the relative positions of waves 1 and 2 at t = 0 resemble those in Fig. 2. Positive values of ω 2 correspond to clockwise rotation. Note that Eq. (11) can represent a wave that propagates in the azimuthal direction λ if w 2 is different than zero. Figure 5 shows the streamfunctions Eqs. (9), (10) and (11) in the horizontal plane for the particular set of parameters indicated in the corresponding caption.
In the panels of Fig. 5 and the following, the centre represents the South Pole and the circular dashed line indicates the Equator. The similarity between Figs. 2 and 5 for the selected set of parameters is evident, taking into consideration that they correspond to stereographic and orthographic projections, respectively. The velocity of fluid parcels in the Cartesian coordinates (x, y) is given by Hamilton's equations: We take the amplitudes to be time dependent in order to emulate changes in magnitudes. Let us start with ε 0 constant and ε 1 = η 1 (1 + sin(µt + π )), ε 2 = η 2 (1 + sin(µt)).
Here, η 1 and η 2 are constants. The time dependence of ε 1 and ε 2 allows us to analyse each wave either separately or together and their transient effect on the observed Lagrangian structures and therefore their transport implications. The time dependence in Eq. (13) is such that one amplitude decreases when the other increases, roughly allowing conservation of the total energy when both waves are present. In the simulations presented below, µ = 2π/10. We begin by considering the case of a mean flow with a = 2 and just wave 2 rotating at different speeds. Furthermore, d = 1 and η 2 = 1. Let us start with ω 2 = 0, i.e. the stationary case. For this case, the dotted line in Fig. 6a shows the azimuthal velocity of the mean flow for ε 0 = −2.5, the dashed line is the azimuthal velocity of wave 2 at λ = 0, where the radial velocity cancels, the solid line is the total azimuthal velocity and the blue line is the wave phase speed. According to Fig. 6a, there are two points where the total velocity cancels, one being the origin. We can also easily see that there are additional fixed points at the r coordinate where the dotted and dashed curves intersect, but placed along the lines λ = π/2, 3π/4. This gives a total of five points in the  hemisphere. Figure 6b shows the M function for τ = 15 evaluated on this stationary field at t = 0. The minima of M highlighting the five fixed points are evident. Moreover, we can see two hyperbolic points in the outer part of the vortex. Next, we consider the case with the same parameters except for ω 2 . Figure 6c shows how this picture changes when ω 2 = 0.1, i.e. for the slow rotation rate of wave 2. The total azimuthal velocity of the wave, in this case, is given by the dashed line in Fig. 6a plus the phase velocity represented by the green line in the figure. If this total azimuthal velocity of the wave is added up to the mean flow, two points are found in which the total azimuthal velocity cancels. Additionally, for a slow rotating wave, similarly to the previous case, the total azimuthal velocity of the wave can still be equal to the zonal mean velocity at some points in the domain. Therefore, Fig. 6c is similar to Fig. 6b except for a rotation. However, for a fast rotation of wave 2 (ω 2 = 4π ; red line), the total azimuthal velocity of the wave will be larger than the zonal mean velocity at all points in the domain. In this case, the pattern of M (Fig. 6d) is very different from the pattern in Fig. 6b showing no HTs. Figure 7 displays the function M obtained from the kinematic model for the same mean flow of Fig. 6a and different parameters for waves 1 and 2. All the representations are for t = 0 and τ = 15. Figure 7a is for the same case as Fig. 6b, except that the amplitude of wave 2 changes in time (η 1 = 0, η 2 = 1). Again, two HTs are visible in the external jet boundary along which filamentation occurs. Figure 7b corresponds to just wave number 1, changing amplitude in time (η 1 = 1, η 2 = 0). We can see one HT at the outer boundary of the jet where material of the vortex is being ejected. In these figures, transport processes producing filamentation-ejecting material have close connections to those present in Figs. 1 and 4a), which have been linked to Rossby wave breaking at midlatitudes (Guha et al., 2016). In Fig. 7c, the mean flow is perturbed by the not-rotating wave 2 of Fig. 7a and wave 1 of Fig. 7b (η 1 = 1,η 2 = 1). In Fig. 7d, the parameters are the same as Fig. 7c, except that wave 2 rotates (ω 2 = 2π/15). The jet shape and filamentary struc- tures greatly resemble those present in the reanalysis data, as shown in Figs. 1 and 4a. Figure 8 present a jet which in the interior is eroded by waves 2 and 1, respectively. To achieve such a configuration, free parameters are specifically tuned, including a zonal mean flow with negative velocities near the pole. In Fig. 8a, the mean flow obtained with parameters ε 0 = 2.6 and a = 0.75 is perturbed by just a travelling wave 2 (η 1 = 0, η 2 = 1, ω 2 = −4π/25) with d = 2. Two filaments projecting material from the interior of the vortex are observed, and they are related to the presence of interior HTs. In Fig. 8b, the mean flow is obtained with the parameters ε 0 = 2.5 and a = 0.5. This mean flow is perturbed by just wave 1 with amplitude that varies in time (η 1 = 1, η 2 = 0). A protruding material filament from the interior of the vortex is observed, which is related to the presence of an interior HT. The interior filaments in these figures recover features that are identified as interior Rossby wave breaking phenomena in de la Cámara et al. (2013) and Guha et al. (2016), and are also visible from the reanalysis data, as shown in Figs. 1 and 4a. Figure 4b shows the pinching of the SPV in the observations on 24 September 2002, which is before its breakup. In the kinematic model, this structure can be obtained with a strong 2 and a substantial contribution from 1 to have a displacement from the pole. Movies S1, S2, S3 and S4 in the Supplement illustrate such structures. In order to reproduce the splitting, we do not need to consider the displacement, and thus we neglect mode 1 in what follows. Figure 9 shows a sequence of M patterns obtained with the amplitude of mean flow is given by where η 0 = −2.5 and µ = 2π/10, and a stationary wave 2 (ω 2 = 0) with amplitude given by Eq. (13). Note that, in this way, the mean flow weakens as wave 2 strengthens, and vice versa. The parameters fit a streamfunction which at t = 0 coincides with that used in Fig. 7a. The development of an hyperbolic point at the pole in the observations (Fig. 4b) can be clearly seen in Fig. 9a. The two vortices have completely split at t = 6.

Kinematic models and conservation of potential vorticity
In this section, we discuss the connection between the kinematic model introduced in the previous section and a fundamental dynamical principle of geophysical fluids. Geophysical flows that are adiabatic and frictionless conserve the po-  tential vorticity Q along trajectories. Conservation of Q is expressed as follows: Here, d /dt stands for the material derivative. A natural question here is to discuss whether the proposed kinematic model conserves Q. Let us assume that our setting is described by the quasi-geostrophic motion of simple vortices in a shallowwater system (see Polvani and Plumb, 1992; Nakamura and Plumb, 1994) in which Q is given by Here, f 0 is a constant related to the rotation rate, D is the mean depth of the shallow-water system, D − h is the total depth, h is the bottom boundary of the fluid layer, which is small when compared to D, and γ = f 0 / √ g 0 D, where g 0 is the gravity constant. is the geostrophic streamfunction for the horizontal velocity field, in our case given by Eq. (7), with parameters corresponding to those of Fig. 7d, i.e. ε 0 = −2.5, η 1 = 1, η 2 = 1, a = 2, d = 1 and ω 2 = 2π/15.
We assume that at the initial time, t = 0, the vorticity Q consists of a circular patch with constant vorticity Q 0 inside and vorticity Q 1 outside. At a later time, t = 2, the vorticity distribution that preserves Eq. (15) is obtained by advecting the circular contour at t = 0, according to motion Eq. (12), with algorithms described in Mancho et al. (2004). Figure 10 summarizes the evolution of the vorticity.
In order to preserve Eq. (16) from time t = 0 to time t = 2, and assuming the barotropic approach in which γ = 0, h is solved from Eq. (16) as Q 0 = 2, Q 1 = 1.8 and f 0 = 20. We note that this calculation could have been repeated for any initial distribution of Q defined as a piecewise constant function. The lower boundary h is thus a time-dependent function adjusted to preserve the conservation of the potential vorticity. Without this forcing, kinematic models would not preserve potential vorticity.

Conclusions
In this work, we propose a simple kinematic model for studying transport phenomena in the Antarctic polar vortex. We are interested in gaining insights into the processes which carry material outward from the vortex structure and inward to the vortex structure. The construction of the kinematic model is realized by analysing geopotential height data produced by the ECMWF. In particular, our focus is on the stratospheric sudden warming event that took place in 2002, producing the pinching and then breaking of the stratospheric polar vortex. We identify the prevalent Fourier components during this period, which consist of a mean axisymmetric flow and waves with wave numbers 1 and 2. The kinematic model is based on analytical expressions which recover the spatial structures of these representative Fourier components. The model can be controlled so that waves with wave numbers 1 and 2 can be switched on and off independently. We are also able to adjust the relative position of the waves with respect to the mean axisymmetric flow.
The study of Lagrangian transport phenomena in the ERA-Interim reanalysis data by means of Lagrangian descriptors highlights hyperbolic trajectories. These trajectories are Lagrangian objects "seeding" the observed filamentation phenomena. The Lagrangian study of the kinematic model sheds light on the role played by waves in this regard. The model is adjusted to a stationary case which considers a mean flow and a stationary wave 2 that perturbs the mean flow in its outer part, producing hyperbolic trajectories. For the stationary case, hyperbolic trajectories are easily identified. This framework is modified by transforming it to a time-dependent problem by making the wave phase speed different from zero or by introducing time-dependent amplitudes. This allows to relate the time-dependent structures with those easily identified in the stationary case. The setting is repeated with wave 1 and both wave 1 and wave 2 together. The joint presence of these waves produces complex Lagrangian patterns remarkably similar to those observed from the analysis of the complex reanalysis data and confirms the findings discussed by Guha et al. (2016). Further adjustment of some model parameters is able to produce erosion by means of filaments just in the interior part of the flow. Finally, we point out that our analysis shows that the breaking and splitting of the polar vortex is justified in our model by the sudden growth of wave 2 and the decay of the axisymmetric flow.
Data availability. The data used in this work are described in Sect. 2.1, where links are also provided to the official websites from which they have been downloaded.