A Three Dimensional Lagrangian Analysis of the Smoke Plume From the 2019/2020 Australian Wildfire Event

During the 2019/2020 Australian bushfire season, intense wildfires generated a rising plume with a record concentration of smoke in the lower stratosphere. Motivated by this event, we use the atmospheric wind reanalysis model ERA5 to characterize the three dimensional atmospheric transport in the general region of the plume following a dynamical system approach in the Lagrangian framework. Aided by the Finite Time Lyapunov Exponent tool (FTLE), we identify Lagrangian Coherent Structures (LCS) which simplify the three‐dimensional transport description. Different reduced FTLE formulations are compared to study the impact of the vertical velocity and the vertical shear on the movement of the plume. We then consider in detail some of the uncovered LCS that are directly relevant for the evolution of the plume, as well as other LCS that are less relevant for the plume but have interesting geometries, and we show the presence of 3D lobe dynamics at play. Also, we unveil the qualitatively different dynamical fates of the smoke parcels trajectories depending on the region in which they originated. One feature that had a pronounced influence on the evolution of the smoke plume is a synoptic‐scale anticyclone that was formed near the same time as, and close to the region of, intense wildfires. We analyze this anticyclone in detail, including its formation, the entrainment of the smoke plume, and how it maintained coherence for a long time. Transport paths obtained with the inclusion of the buoyancy effects are compared with those obtained considering only the reanalysis velocity.

• The impact of the vertical velocity and the vertical shear is described using different Finite Time Lyapunov Exponent (FTLE) formulations • Coherent regions where smoke parcels have qualitatively different fates are characterized with and without considering buoyancy effects • Lagrangian Coherent Structures linked to the anticyclone, which affected the smoke plume's evolution, are analyzed in detail Supporting Information: Supporting Information may be found in the online version of this article.
The goal of this paper is twofold: first, we are interested in applying the LCS approach to better understand some aspects of the atmospheric circulation in general, and second, we want to explore the link between some of the uncovered LCS and the smoke plume evolution.We divided our study on three periods of time representative of the three phases of the smoke plume event: the start of the first pyroCb event (late December); the time when the cloud reaches a highest concentration of smoke (early January), and almost 2 months later to see how the particles have persisted in the stratosphere (late February-early March).

Data
Our Lagrangian analysis is based on the ERA5 reanalysis data set, the fifth generation of the European Center for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalysis produced by the global climate Copernicus Climate Change Service (C3C) (Hersbach et al., 2019).The ERA5 provides lateral wind velocity    (m/s), vertical velocity ω (Pa/s), geopotential, and temperature in 37 pressure levels from 1,000 to 1 hPa.The temporal resolution is 1 hr.
The ECMWF provides the vertical velocity ω in Pa/s with negative values corresponding to upward motion.To compute the vertical velocity in meters per second we use the hydrostatic approximation which assumes that the horizontal scale is larger compared to the vertical one, that is, where ρ is the density, g is the gravity, and w is the vertical velocity in m/s.Here, the density is related to pressure p and the temperature T through the equation of state of the ideal gases, p = RρT with R = 287.058(m 2 /s 2 K −1 ).
In order to track the movement of the real observed plumes of smoke in the stratosphere we use the Earth Observing System (EOS) Aura Microwave Limb Sounder (MLS) Level 2 standard product for geopotential height.The data version used is 5.0 (Schwartz et al., 2020).MLS provides day and night near-global (82°S-82°N) measurement of vertical profiles of various atmospheric gaseous compounds geopotential height, and temperature of the atmosphere.The measurements yield around 3,500 profiles per day for each species with a vertical resolution of approximate 3-6 km.Following Kablick et al. (2020), we use the information of the water vapor mixing ratio H 2 O, the collocated carbon monoxide mixing ratio CO and the geopotential height from the MLS data set.The

Lagrangian Methods
We work in the Lagrangian framework, that is, we analyze stratospheric transport following parcels' trajectories.
We look for LCS which control the stretching and folding of the polluted air mass and separate regions where trajectories have qualitatively different kinematic fates.The method that we use to approximate LCS is the FTLE (Haller, 2000(Haller, , 2002;;Shadden et al., 2005) which measures the exponential separation rate between initially nearby air parcels.
Let x(t; x 0 ) be a trajectory of an air parcel that starts at x 0 at time t 0 : where v(x, t) is the velocity vector field.Thus,

𝐯𝐯𝑑𝑑𝑡𝑡𝑑
Let F be the strain tensor given by and G = F ⊺ F be the right Cauchy-Green tensor.The FTLE is defined by where λ max (G) is the maximum eigenvalue of the matrix G. Repelling LCS are defined as maximizing ridges of the FTLE field computed from forward trajectories (final time t f > t 0 ), and attracting LCS are defined as ridges in the backward-time (t f < t 0 ) FTLE field.The repelling and attracting LCS identified in this manner are proxies for the finite-time counterparts of the stable and unstable manifolds of hyperbolic trajectories from the classical theory of dynamical systems (Balasuriya et al., 2018;Branicki & Wiggins, 2010;Haller, 2011).Thus, FTLE field are usually considered as indicators of hyperbolic LCSs (Haller, 2001), although they can produce both false positives, where separation is due to shear and not hyperbolic behavior, and negatives in LCS detection (Haller, 2002(Haller, , 2011) ) even in simple two-dimensional steady flows (Farazmand & Haller, 2012;Haller, 2011).Note, however, that because linear shear yields linear separation between particles and linear strain yields exponential separation, in most realistic oceanic and atmospheric flows, where hyperbolic strain-dominated regions are in abundance, separation between trajectories is indeed dominated by hyperbolic behavior, and so FTLE ridges very rarely produce false positives.In our case FTLEs seem to be a useful tool that produces physically relevant partitions of the domain.This is further tested and confirmed for all the FTLE ridges that we present in this paper by considering evolution of trajectories originating on opposite sides of the FTLE ridges.For all of our ridges, trajectories on opposite sides exhibit qualitatively different fate, and the separation between them is not aligned with the FTLE ridge (as it would have been for shear-dominated ridges).

Exact and Reduced Formulations of FTLEs, and Their Use as a Quantitative Measure of the Relative Importance of Vertical Velocity and Vertical Shear in the Atmospheric Flows
In geophysical flows the full computation of the evolving 3D velocity field is challenging, and the vertical velocity, w, which is generally much smaller than the horizontal velocities, is often estimated as a diagnostic quantity 10.1029/2023JD039773 4 of 21 (rather than prognostically solved as part of the equations of motion like the horizontal velocity components) (Donald Ahrens & Henson, 2015;Hersbach et al., 2019) and is thus less reliable.It is hence tempting to ignore w in the computation of FTLEs.However, as we show below, this approach leads to large errors in situations where vertical shear of horizontal velocity is large.
In order to investigate the effects of the vertical velocity and vertical shear on the resulting FTLEs and LCS, following Sulman et al. (2013), we compare the reduced FTLE formulations given by the definitions: Case 1. 2D form of tensor G with trajectories from 2D, Case 2. Case 1 using 3D trajectories   = (  ) ∈ ℝ 3 , that is, Case 3. 2D form of tensor G with vertical velocity and trajectories from 3D, Case 4. 2D form of tensor G with vertical shear and trajectories from 2D, that is, Case 5. Case 4 using 3D trajectories   = (  ) ∈ ℝ 3 , that is, To quantifying the effects of vertical velocity and vertical shear on the resulting spread of trajectories (i.e., FTLEs) we perform several trajectory calculations.First, we advect trajectories using only the horizontal components of velocity, that is, advection along a constant height surface, and we refer to these trajectories as 2D.Next, we advect trajectories starting at the same initial conditions using all three components of velocity, and we refer to those as 3D.Please note that the temporal evolution of trajectories in 3D modifies the height according to the 5 of 21 corresponding time and space dependent w − component of velocity.Then, we compare the different formulations of FTLE using air parcel trajectories restricted to 2D (Case 1 and 3) or allowing the air particles to move in three-dimensional space (Case 2) and also with and without the terms in F associated with the vertical shear (third column) and/or vertical velocity (third row).
For the computation of the FTLE, trajectories are determined by integrating the differential Equation 2 using a fifth-order Runge-Kutta method (Cash-Karp) with a fixed step size of 1 hr, which provides estimates accurate to fifth order.The derivatives in the Cauchy-Green tensors are then approximated using second-order centered finite-differences and the eigenvalues are calculated with the MATLAB function eig that use the QZ algorithm also known as the generalized Schur decomposition.For the FTLE study in this section, trajectories are estimated over a time interval of 5 days, which is sufficiently long for the ridges in FTLE fields to become well-defined, but sufficiently short to not produce overly complex and tangled ridges.
The stable manifolds (repelling structures) are calculated through FTLE using forward trajectories.The unstable manifolds (attracting structures) are calculated through FTLE using backward trajectories.In the figures below, we represent repelling structures in blue and attracting structures in red.
Since we are working in the atmosphere, we change our coordinate system from Cartesian to spherical.Therefore the matrix F has the following form, where r is from the Earth center height, θ is a function of the latitude and φ is longitude in radians.
Section 3 is devoted to quantifying the effects of vertical velocity and vertical shear on the resulting spread of trajectories (i.e., FTLEs).For this purpose, we perform several trajectory calculations.First, we advect trajectories using only the horizontal components of velocity, that is, advection along a constant height surface, and we refer to these trajectories as 2D.Then we advect trajectories starting at the same initial conditions using all three components of velocity, and we refer to those as 3D.Please note that the temporal evolution of trajectories in 3D modifies the height according to the corresponding time and space dependent w-component of velocity.In the rest of the manuscript, all trajectories are advected using all three components of velocity, that is, they move in 3D.
Figure 1 shows the different cases of forward FTLE (first column) and backward FTLE (second column) computed with τ = 5 days at 11 km height.The initial condition for this calculation is 22 December 2019, with an initial height of 11 km.Case 1 and Case 2 are too low in magnitude compared to Case 6, indicating stronger lateral separation of trajectories caused by the 3D flow over the same time interval.Note that, although magnitude of FTLEs is important as a measure of separation, it is the location of the FTLE ridges, not their strength, that is relevant for identifying LCSs.Comparing the geometries of FTLE regions, we observe that cases 1 and 2, which ignore w and vertical shear, produce similar but slightly displaced large-scale LCS compared to Case 6, but grossly under-estimate small scales.The main differences between cases 1-2, which ignore the vertical velocity and the vertical shear, and six are underestimate small scales and the large-scale (area-averaged) offset (i.e., the difference in magnitude).This seems to be less important in the subtropical jet region south of about 40°S, which is dominated by the large-scale LCS, but it seems to lead to large discrepancies north of 40°S, including region over Australia, where small scales nearly erase the large-scale ridges seen in the upper two rows of Figure 1.Case 3, which considers vertical velocity but ignores vertical gradients, is also too low in magnitude showing the same issues as Case 2 in terms of producing overly strong and slightly misplaced large-scale FTLE ridges and missing nearly all small-scale ridges.On the contrary, Case 4 which includes the effects of vertical shear in the horizontal velocity components, improves in magnitude, and better represents the larger-scale LCSs in terms of both their strength and location, and starts to capture some (but not all) small scale features.However, many of the small scales over Australia and in the northern part of the domain are clearly still missing.Case 5 is pretty close to Case 6 (i.e., full 3D formulation of FTLEs) in terms of both magnitude and ridges locations.Note that, the difference of Case 5 with Case 6 is five orders of magnitude smaller than the differences for the other cases.The drastic improvement of Case 5 compared to Case 3 suggests that the influence of vertical shear on the spread of 3D trajectories is much more important than the influence of vertical velocity.The significant improvement of Case 5 over Case 4 highlights the significant differences between the lateral spread of 3D trajectories versus 2D trajectories and thus points to the importance of using 3D trajectories in the computation of FTLEs.It is thus extremely important to both use 3D trajectories and include the terms corresponding to the vertical shear in the computation of FTLEs.This is similar to the situation in the ocean at submesoscale, but in ocean mesoscale flows, Case 4 is typically closer to Case 5 (Lanotte et al., 2016;McWilliams, 2016;Sulman et al., 2013).The same cases are studied for 5 km height (lower troposphere) and 20 km (stratosphere) in the Supporting Information S1 (see Figures S1 and S2).Similar conclusions are reached.
In the upper troposphere and stratosphere the vertical velocity in isentropic coordinates is small, except near active convection.Therefore, a widely used method to calculate trajectories in the atmosphere is to assume that the particles are constrained to remain on surfaces of constant potential temperature (θ). Figure S4 in Supporting Information S1 shows that computing 3D trajectories using three components of velocity in z-coordinates is consistent with advecting trajectories along isentropic surfaces.Therefore, either of these methods could be used to compute air parcel trajectories and estimate FTLEs at any given location.In the remainder of this paper, we use the first method.A proper rendition of the 3D FTLE fields can then be achieved by either vertically stacking together FTLEs at different z-levels or stacking together FTLEs at different isentropic levels.In the remainder of this paper, we will do the former.

Geometry of the LCSs in the Stratosphere Near Australia During the Aerosol Injection
The first relatively small pyroCb event occurred around 22 December 2019 (Peterson et al., 2021), and the main smoke ejection event occurred on December 31.Using the OMPS AI from NASA's Fire Information for Resource Management System (Flynn et al., 2014), panel (a) of Figure 2 shows the Absorbing Aerosol Index highlighting the smoke plume generated by this first aerosol injection.Similarly, panel (b) shows the plume generated by the In this section we mapped out LCS over and around Australia on late December.Our goal here was two-fold: first, we wanted to see what sort of the 3D LCS geometries existed in the stratosphere at the time of first smoke injection; and second, whether any of these structures were influencing the movement of the actual smoke plume.The first question is more of a generic study of possible 3D LCS geometries in the stratosphere, while the second question is more applied.
Starting with the first question, we show in Figure 3 the horizontal slice of 3D forward FTLEs on December 22 at 11 km height.Figure 3a shows where the main aerosol plume represented in Figure 2a on December 22 is located with respect to the LCSs.Multiple FTLE ridges can be identified over and near Australia, and we will next consider three of them (highlighted by the black, red and blue arrows in Figure 3b) in some detail.We specifically picked these three ridges because they have qualitatively different geometry in 3D.
The black arrow in Figure 3b points to an elongated nearly zonal ridge slightly south of 40°S that is folded onto itself near 120°E.It is located near the northern edge of the subtropical jet.The red arrow points to another nearly zonal ridge just north of 40°S, which cuts off the very southernmost tip of Australia near 145°E.And the two blue arrows indicate what seems to be either one continuous ridge that is folded near 150°E, or perhaps two distinct ridges that come close to each other near 150°E.
Although in the horizontal slice (Figure 3b) all three FTLE ridges seem to be rather similar (all look like wiggly 1D curves in the horizontal slice), they look remarkably different in the vertical slice (Figure 3c shows a vertical slice of 3D FTLEs at 140°E).The same three features are marked by same-colored arrows.
The simplest geometric structure (see red arrows on panels (b and c)) resembles a vertically tilted curtain spanning about 8 km in vertical, from ∼5-13 km.In order to gain some insight into which qualitatively different regions this barrier separates, we released trajectories on either side of it (marked by red dot to the north and red square to the south in Figures 3b and 3e).South of this barrier, parcels move rapidly to the east, generally maintaining or even increasing their altitude (with a bit of altitude decrease at the end of 5 days).North of this barrier, the parcels also move to the east, but with more northward deviation, much more slowly, and at a generally lower altitude (see an example of the trajectories in Figure 3e; colors represent the height of the trajectory in km).
The LCS pointed by the blue arrows in Figures 3b and 3c has a slightly more complex vertical structure.This structure, shaped like a hat, acts as a lid preventing upwards vertical transport.This cap-like structure divides Australia into three regions as shown in panel (b) and (c).Parcels that originate outside the hat (i.e., to the north of the northern segment of the ridge and south of the southern segment ridge highlighted in Figure 3b or above it) move eastward increasing or maintaining height.However, parcels that originate inside/underneath the hat structure move west and down, as shown in 3D.
Finally, the LCS marked with black arrow has the most interesting geometry out of the three.Topologically, it is a tube (and thus appears as a closed curve in a vertical slice in Figure 3c i.e., closed at its western end where the manifold in Figure 3b folds onto itself creating a closed elbow.An example of two trajectories inside and outside the tube is shown in Figure 3f but the geometry in the neighborhood of this tube-shaped structure is more complex and requires a more in-depth study.As we will see in the next figure, this tube structure seems to also interact with the nearby eddy located to the northeast of it, and with another eddy located further to the southeast. First and second column of Figure 4 show four daily snapshots of the forward (blue) and backward (red) FTLEs at 11 km height near the tubular structure from 22 to 25 of December.The third column of Figure 4 is a schematic diagram showing the intersecting attracting and repelling LCS near the tube.The Lagrangian geometry in this region is governed by two hyperbolic trajectories (HT1 and HT2) that give rise to two pairs of intersecting stable and unstable manifolds.(The stable and unstable manifolds of HT1/HT2 are shown in purple/blue and red/orange).The tube is nothing other than a lobe that is trapped by a segment of unstable manifold of HT1 and a stable manifold of HT2.Initially, this lobe is close to HT1 but moves toward HT2 with time.As it does so, the segment of its bounding unstable manifold elongates, and the segment of the stable manifold shrinks, so as the tube moves away from HT1, it gets shorter and wider.Later on, as it approaches HT2, it becomes stretched along the unstable manifold of HT2 and becomes narrow and long again.This is a classical picture of a heteroclinic tangle, which suggests that the turnstile lobe mechanism is a common phenomenon in the stratosphere.The presence of lobe mechanics in the stratosphere has been showed by several authors (Joseph & Legras, 2002;Koh & Legras, 2002;Koh & Plumb, 2000), although in those previous cases the lobes involved in the lobe turnstile mechanism had a curtain-like geometry (more akin our red FTLE ridge in Figuresure 3b and 3c), rather than tubular geometry.
The behavior of the different sets of particles in and around the tube-lobe is also shown in Figure 4.The black and green particles are released inside and outside the tube, respectively.The purple dots correspond to particles that Green color identifies parcels that on December 22 are outside the tubular structure formed by the stable manifold but close to it.Purple color identifies parcels on the eddy structure and black color identifies parcels inside the tube.The third column shows a diagram of the relative position of the stable (in purple and blue) and unstable (in red and orange) manifolds associated with the hyperbolic trajectory HT1 and HT2, respectively.are released inside the nearby northeastern eddy.Black particles move eastward following the stable and unstable manifold of HT1. Green parcels also move eastward but following the stable manifold of HT2.On 24 Dec 2019, both set of particles approach HT2, and their route is interrupted by the unstable manifold of the HT2 (orange line in the diagram).From there on, black and green parcels diverge and move in different directions, as seen in Figure 4 (December 24-25).
Regarding the second question, the presence of the three structures studied in this section did not play a crucial role in the smoke plume's evolution, because the actual smoke plume on December 22 mainly falls on top of the area with many tangled FTLEs ridges east of Australia, as shown in panel (a) of Figure 3.Such regions mark areas of rapid stirring and mixing which suggests that the plume overlaying that region on December 22nd will disperse and will be unlikely to stay coherent for long.To provide a complete analysis of late December, we have also included a study of the LCS in relation to the position of the smoke plume on December 31 (the day of the main injection) in the Supporting Information S1 (see Figure S5).While the LCS described in this section did not aid in explaining the smoke plume's later evolution, they contribute to a more comprehensive understanding of the atmosphere's structure and geometry before the event.Features such as an anticyclone that influenced the evolution of the smoke plume were actually produced a few days after the main injection.The LCS associated with these features will be described in detail in the next section.

Split of the Main Aerosol Plume (Early January 2020)
The observed plume was visible from satellites starting shortly after the injection, although with limited resolution.Over 1-4 January parts of the plume were detected moving to the southeast from Australia, and on January 6th a very coherent patch was identified near 120°W; 50°N (Kablick et al., 2020).
Figure 5a shows the observed plume, as detected from the satellite, on January 6th and onward to end of February.On January 6th, the highly concentrated plume reaching roughly 1,000 km in diameter was detected in the stratosphere at about 15 km near 100°W, 60°S.From there, the smoke plume split into three parts that moved along three different paths (Kablick et al., 2020).Path P3 (blue dashed line in Figure 5a) went eastward at a nearly constant height of about 17 km, whereas path P1 (green solid curve in Figure 5a) looped around and went westward, ascending on its way and passing south.P2 drifted toward the south over Antarctica on January 5 (red dashed line) and remained south of 60°S until mid-February.Path P2 not studied here.
Motivated by the split of the plume on January 6th and the striking difference between the P1 and P3 paths, we have decided to focus in P1 and P3 and applied the Lagrangian approach to better understand the cause of this splitting, the subsequent transport geometry, and the influence of the plume buoyancy on its movement.
Consistent with observations, simulated parcels in the ERA5 model released in the area of the observed plume (black box in Figure 5) on Jan 6th also showed the splitting into two distinct P1-like and P3-like groups.Without the buoyancy effects, however, the P1 path is too low in altitude, is shifted northward, and passes over Australia rather than south of Australia on February 26, as in observations.Using the difference in altitude between the simulated and observed P1 trajectories on February 26th, we have estimated the time-averaged buoyant velocity to be about 0.0022 m/s.When this buoyant velocity was added to the ERA5 velocities, the agreement with observations significantly improved.With buoyancy, P1-like path shifted up and south, with P1-like trajectories passing south of Australia by February 26th, consistent with observations.Advecting P1-like trajectories backward in time (with negative buoyancy of 0.0052 m/s) from January 6, we observed that these passed just to the southeast Australia on December 31st, that is, right within the area of the observed plume of the main ANY event (see Figure 2b right).This suggests that the plume observed on January 6th was likely generated by the main event on December 31, rather than the earlier event of December 22.This also agrees with our previous conclusion that the plume generated on December 22nd was unlikely to stay coherent for long.
To simplify the analysis, we assumed a constant buoyancy during the ascent of the smoke.However, it is important to note that treating buoyancy as the sole driving force for upward movement represents an oversimplification of the complex physics involved in plume rise.This approach could mask the underlying complexities and may not capture the nuances of the actual processes (Davies, 2015;Khaykin et al., 2020;Lestrelin et al., 2021).
The splitting of the plume into P1 and P3 on January 6 suggested the presence of strong LCS in the area at that time, which acted as transport barriers with different trajectory fates for parcels on the opposite sides of LCS.CURBELO AND RYPINA 10.1029/2023JD039773 14 of 21 centered at around 125°W, 55°N with two long and narrow tendrils extending from it.The latter originate in the eastern part of the domain and are separated from the rest by a strong FTLE ridge, which has a tilted-curtain-line geometry in 3D spanning the altitudes of the observed plume (16-22 km).Panel (d) of Figure 6 displays two different angles of a tilted curtain that separates parcels which move westward (left and above the blue structure on January 6) or eastward (right and below the structure).LCS geometry is qualitatively similar to that shown in Figure 6 at altitudes ranging from 16 to 22 km.
It is interesting to look at the changes in LCS geometry with the addition of buoyancy.Without buoyancy, the eddy only contains a small percentage of P1-like trajectories (red), with most parcels (yellow) continuing westward, rather than eastward after looping around.With buoyancy, however, almost the entire eddy becomes red.This is because buoyant parcels rise up higher and are then carried westward by the strong westward winds at higher altitudes.
In the previous studies, there are two views on the origins of P1 and P3 trajectories.(Kablick et al., 2020) suggests that the main smoke plume has split in P1 and P3 on 6 January 2020 (Khaykin et al., 2020), however, argues that P1 and P3 have different origins; P1 was produced by the smoke ejection on the 31st December, whereas P3 was originating from a different smoke ejection event on January 4th.Our analysis show that on January 6, in the general area where smoke plume was detected, all P1-like trajectories were located inside an anticyclonic eddy.Additionally, all P3-like trajectories were located about 30° further east from the center of the eddy, on the eastern side of a separating manifold (see Figure 6).Furthermore, our analysis indicates that when these trajectories were advected backward in time from January 6 to January 4th, and then to December 31, both P-1 and P-3-like trajectories took similar paths and backtracked to the same general area south-southeast of Australia (see Figure 5d).The accuracy of our Lagrangian calculations does not definitively determine whether the P-1 and P-3 trajectories were generated by the same smoke event on the 31st or by two distinct smoke injection, one on the 31st and another a few days later.

The Anticyclonic Vortex Evolution
The anticyclonic vortex, which contained all P1-like trajectories on January 6 in our simulations, has been detected and associated with the smoke bubble using observational data in a recent paper by Khaykin et al. (2020).Specifically, in that paper, the plume started to encapsulate January 4 and the coherent smoke bubble was first observed on January 7 and was traced continuously to early April.
Motivated by this remarkable observation, in our subsequent analysis we use LCS techniques to trace the anticyclone evolution from its formation around January 5th to its eventual breakdown in mid-March.When mapped against available observational data, we will see that a coherent bubble of the smoke was always confined by the eddy throughout its entire lifetime.Figure 7 is intended to illustrate various aspects related to the eddy, including its formation, evolution, entrain of the smoke plume, coherence maintenance of the anticyclone, and confinement of the smoke within it.In what follows, we also try to shed light on the shedding of material from the anticyclone during its ascent into the stratosphere and the eventual breakdown of the anticyclone.
In Figure 7a, we map out LCSs starting on 29 December through 6 January.As this figure illustrates, close to its formation near January 5th, the anticyclone was initially part of the dipole structure.Specifics of its formation and initial evolution are shown in inset panels of Figure 7a and can be described as follows.On December 29th, a cyclonic eddy was shed by the stratospheric jet in the region southeast of Australia near 180°E and 55°S, not too far (∼1,000 km) from where the smoke plume was observed on December 31st.The cyclonic eddy started propagating to the northeast, reaching 190°E by December 31 and 240°E by 4 of January.Around the same time, an anticyclonic counterpart has started forming just south from the cyclonic eddy, with the classic dipole geometry fully formed, that is, closely bound pair of vortices of equal strength with opposite circulation, and clearly visible on 5 and 6 of January.Shortly after that, the cyclonic eddy started to disintegrate quickly and was completely gone a few days later.The anticyclonic eddy, however, remained.This anticyclonic part of the dipole is indeed the same eddy that has been shown in Figure 6 to contain P1-like trajectories, and which we will show a confined smoke bubble for the next two-three months (Figures 7b and 8).
Overlaying available satellite observations of H 2 O vapor (indicator of smoke plume, black and empty circles in Figure 7a) on top of LCS structures reveals that the presence of smoke plume in and around the anticyclonic eddy on January 5.It is this smoke that got entrained into this eddy shortly after its formation that will later take the P1-path in Figure 5 and will largely remain in the eddy until the beginning of March.
In Figure 7b, we now follow the anticyclonic eddy throughout its lifetime from 6 January to 4 March (the time it starts to decay), mapping the eddy-defining LCSs against available ozone anomalies observations (another CURBELO AND RYPINA Was the smoke confined to the eddy in the vertical as well as in horizontal?The answer is yes, as seen from mapping out maximum ozone anomalies concentration on top of LCS slices at different heights (ozone anomaly is from Khaykin et al. (2020)).Specifically, in Figure 8 we clearly observe a great correspondence between ozone at a given height (empty red dot means no or very low ozone anomaly and filled red dot means high ozone concentration) and the eddy LCS signature (filled red dot is always located at the center of strong FTLE ridges delineating the periphery of the eddy).
On 26 February the smoke plume reached above 31 km (see panel (a) in Figure 5, path P1).Regarding our simulated plume, the locations of the forward-tracked (from January 6 to February 26) air parcels (red dots) over the backward FTLEs (proxy attracting LCS) on that day are displayed in Figure 9a.Both FTLEs and trajectories were computed here including buoyancy effects.Although most smoke parcels remained within the eddy, some smoke parcels were slowly escaping, or leaking out of the eddy, along the unstable manifolds.The parcels that leaked out were rapidly stretching into long filaments and were then mixed rapidly with surrounding air so that the concentrations became small.Panel (b) of Figure 9b shows the superimposition of the MLS aerosol data on top of the attracting LCS.The formation of the filaments associated with the leaking vortex was discussed by Khaykin et al. (2020) and Lestrelin et al. (2021), who observed the leakage from the vortex bottom or tail using Cloud-Aerosol Lidar with Orthogonal Polarization data.

Conclusions
The stratospheric winds have relatively weak vertical velocities compared to horizontal velocities.It is thus tempting to ignore w and consider the motion of air parcels in 2D.However, as we have shown using different formulations of FTLEs in this paper, such 2D approach is misleading as it does not take into account the fact that even slight vertical movement might expose air parcels to different horizontal advection due to strong vertical shear.Thus, for an accurate representation of 3D transport, it is necessary to consider the movement of trajectories in 3D, and, importantly, include the vertical shear terms in the formulation of the FTLE matrix.On the other hand, due to the smallness of vertical displacement of trajectories, the terms associated with the vertical movement itself can be safely ignored, so the formulation of FTLEs can be reduced from 3 × 3 matrix to the 2 × 3 matrix without much reduction in accuracy.Note that there is an important conceptual distinction between what we refer to as 2D trajectories (i.e., trajectories computed using horizontal velocity at constant height) and quasi 2D trajectories advected along isentropic levels.The latter, in contrast to the former, are nearly equivalent to 3D trajectories.This is because in large-scale atmospheric flows, geometric altitude z is not the optimal vertical coordinate.In the stratosphere in particular, the natural vertical coordinate is potential temperature θ, which naturally separates the across-isentropic thermodynamic effects of the flow from the along-isentropic flow Motivated by the strong Australian wildfire event in 2019/2020, we have applied the Lagrangian approach to study the 3D transport in the stratosphere.The study is based on the ERA5 reanalysis winds and compares simulations of the smoke plume with available observations.Over the last few decades, FTLEs has been shown to provide a useful tool for mapping out transport properties in geophysical flows.Here, we used FTLEs to uncover several distinct geometries of 3D transport in the atmosphere.Among them were curtain-like 2D sheets of transport barriers, a 2D hat-like structure, and a tubular structure, which upon inspection proved to be a turnstile lobe.The "curtains" have been observed in prior work, as well as curtain-like turnstyle lobes (Joseph & Legras, 2002;Koh & Legras, 2002;Koh & Plumb, 2000).However, we are not aware of the appearance of "hat" and "tube" in prior literature.
In the second part of the paper, we have applied our FTLE-based approach to better understand the evolution of the smoke plume from the wildfire event of 2019-2020.We mapped out LCS on January 6th, the day when a very coherent smoke patch was detected from satellites almost half-way across the globe from Australia.This patch was then observed to split into two parts, one moving eastward at low altitude and another looping around and heading westward back toward Australia at much higher altitudes.Analysis of simulated trajectories in the ERA5 reanalysis model suggested that buoyancy of hot smoke strongly affects the movement of smoke plumes, but the separation of the main smoke plume into two different paths occurs both with and without buoyancy.With buoyancy included, simulated plume matched the observed one very well (without buoyancy, it did not).We then used the simulated trajectories with buoyancy to map out LCS on Januay 6th, which clearly delineated regions destined to take two different paths.We found a tilted surface that divide the region in two areas.Trajectories starting anywhere to the east of that tilted curtain behave move eastward at a lower altitude (like P3), whereas trajectories to the west of the curtain move westward.Thus, the switch from P1 to P3 behavior can be achieved by either shifting location of the particle on the January 6th in horizontal (east to west of the tilted curtain) or in vertical (above to below of the tilted curtain).
Our analysis also suggested that the bulk of the smoke bubble was largely contained by an anticyclonic eddy that formed southeast of Australia shortly after the smoke injection event.This is in line with prior work by Khaykin et al. (2020) and Kablick et al. (2020).
It is interesting that our analysis suggested that the anticyclone was initially formed as part of a dipole structure, but became an isolated single eddy after its cyclonic counterpart decayed in the beginning of January.The anticyclone was rising up in the atmosphere throughout its lifetime, and changed its vertical extent from a few kilometers near its formation to 6 km in mid-February and then back to a few kilometers in early March when it decayed.
When mapped against satellite observations of smoke indicators, the eddy as mapped by FTLE ridges was always aligning well in both horizontal and vertical with the smoke bubble.This supports the earlier conclusion of Khaykin et al. (2020) and Kablick et al. (2020) that it was this eddy that largely shaped the evolution of the smoke bubble.Khaykin et al. (2020) and Lestrelin et al. (2021) highlighted that heating caused by the smoke was a key factor in sustaining and stabilizing this anticyclonic vortex.In the beginning of January, we observed the formation of a dipole, with the smoke present only in the anticyclone part.Additionally, we noted the decay of the cyclonic vortex, which did not contain aerosol inside.However, based on our analysis, we cannot definitively conclude whether the robustness of the anticyclone was due to the presence of smoke or other reasons.Superposition of forward-and bachward-FTLE analysis revealed that, while most of the smoke has been well contained by the eddy core, some smoke was slowly leaking from the eddy (along the backward-FTLE ridges, as expected from dynamical systems theory).It is not entirely clear why this leakiness is not observable from satellites, but perhaps the answer is simply that because of the elongated shape of the leaky filaments (and thus increased tracer gradients across the filament), they would mix rapidly with surrounding air to undetectable concentration levels.
Because of its chemical and physical properties, hot smoke interacts and feedbacks on the atmosphere in a complicated way.Correctly accounting for these effects is challenging, which complicates the prediction of the smoke plume evolution.Here, we have suggested a simplified way in which all of these complex effects could be parameterized by a single parameter, a constant additional vertical velocity (0.0022 m/s as estimated from best fit 10.1029/2023JD039773 19 of 21 to available data), which can be added to the atmospheric circulation model to account for the rising of the smoke plume with time.Simulations of the smoke evolution performed with and without additional vertical velocity were implemented, and it was shown that adding constant vertical velocity greatly improved the correspondence between predicted and observed paths of the smoke.We also backtracked trajectories that clearly point to the December 31st event as the source of the patch, in agreement with other previous satellite imagery (Khaykin et al., 2020) and modeling (Yu et al., 2021) studies.
Overall, this work suggests that FTLEs might be a useful tool in understanding and predicting the evolution of a pollutant patch, specific for the 2019/2020 bushfire event studied here, or more generally in other atmospheric applications.

Figure 1 .
Figure 1.Forward (first column) and backward (second column) Finite Time Lyapunov Exponent (FTLE) computed with τ = 5 days for the different formulations at 11 km above the sea-level.Each row corresponds to one of the FTLE formulations described in the text.To quantify the effects of vertical shear and vertical velocity the third column represents the absolute error between the forward FTLE for the corresponding case to the fully 3D FTLE (Case 6).

Figure 2 .
Figure 2. Ozone Mapping and Profiler Suite Aerosol Index (AI) colormap over Australia during the days (a) 22 December 2019 and (b) 31 December 2019.The red dots are fires/hotspots from VIIRS given by thermal anomalies.UV AI values greater than 5 marks the smoke plume.The black dots are the positions on the corresponding day of the backtracked pollutant parcel trajectories starting on January 6 discussed in Section 4.2.1.

Figure 3 .
Figure 3. Panel (a) shows where the main Aerosol plume (Figure 2a) is located with respect to the LCSs on December 22. Panel (b) and (c) show the horizontal and vertical slices of the 3D forward Finite Time Lyapunov Exponent computed the day 22 December 2019 with τ = 5 days at z = 11 km height and a fix longitude 140° E, respectively.The dashed line in panel (c) represents the tropopause.The arrows highlight the three coherent structures described in the text.The markers correspond to the initial position of parcels trajectories 5 days forward that are shown in panels (d)-(f) as examples of dynamics in the different regions.

Figure 4 .
Figure 4. Snapshot of the time evolution of Finite Time Lyapunov Exponent forward (first column) and backward (second column).Panel displays forward parcel trajectories at 11 km that are initialized on 22 December 2018.Green color identifies parcels that on December 22 are outside the tubular structure formed by the stable manifold but close to it.Purple color identifies parcels on the eddy structure and black color identifies parcels inside the tube.The third column shows a diagram of the relative position of the stable (in purple and blue) and unstable (in red and orange) manifolds associated with the hyperbolic trajectory HT1 and HT2, respectively.

Figure 6 .
Figure 6.(a)-(b) Simplified scheme of forward path types initialized on 6 January 2020 in different regions (in colors) at 18 km height.Panel (a) is obtained from the calculations using the original reanalysis rates (without adding the buoyancy effect).Panel (b) considers a constant buoyancy so a value of 0.0022 m/s is added to the vertical velocity.Forward Finite Time Lyapunov Exponent are represented in black and in gray for τ = 20 and τ = 40 days respectively.(c) Scheme of possible paths by color.The red trajectory corresponds to path P1 and the magenta one to path P3 in Figure 5.(d) 3D curtain that separates the trajectories of particles with behavior similar to P1 with those whose patch is P3 (with buoyancy).
on January 15-26 km on February 9 and 27 km on March 4. Second, the vertical extent of the anticyclone, which is about 4 km close to its formation, seems to first increase to 7 km by February 9 and then decrease to just 3 km by March 4.

Figure 9 .
Figure 9. Longitude-latitude section of Backward Finite Time Lyapunov Exponent (color background) computed in 3D with buoyancy with τ = 40 days on 26 of February 2020.In panel (a), the red dots are the final position of the trajectory that started on January 6 represented in the Figure 5(b).In panel (b), dots display the tracking position of the H20 maximum Microwave Limb Sounder profile plume for the corresponding day.Filled blue circles represent values that are inside the 1 km interval centered at the corresponding height level.Empty blue circles represent values that are outside that interval.