Lower Thermospheric Material Transport via Lagrangian Coherent Structures.

We show that inter-model variation due to under-constraint by observations impacts the ability to predict material transport in the lower thermosphere. Lagrangian coherent structures (LCSs), indicating regions of maximal separation (or convergence) in a time-varying flow, are derived in the lower thermosphere from models for several space shuttle water vapor plume events. We find that inter-model differences in thermospheric transport manifest in LCSs in a way that is more stringent than mean wind analyses. LCSs defined using horizontal flow fields from the Specified Dynamics version of the Whole Atmosphere Community Climate Model with thermosphere-ionosphere eXtension (SD-WACCMX) at 109 km altitude are compared to Global Ultraviolet Imager (GUVI) observations of the space shuttle main engine plume. In one case, SD-WACCMX predicts an LCS ridge to produce spreading not found in the observations. LCSs and tracer transport from SD-WACCMX and from data assimilative WACCMX (WACCMX + DART) are compared to each other and to GUVI observations. Differences in the modeled LCSs and tracer positions appear between SD-WACCMX and WACCMX + DART despite the similarity of mean winds. WACCMX + DART produces better tracer transport results for a July 2006 event, but it is unclear which model performs better in terms of LCS ridges. For a February 2010 event, when mean winds differ by up to 50 m/s between the models, differences in LCSs and tracer trajectories are even more severe. Low-pass filtering the winds up to zonal wavenumber 6 reduces but does not eliminate inter-model LCS differences. Inter-model alignment of LCSs improves at a lower 60 km altitude.

Introduction

Material transport in the lower thermosphere is of scientific interest. Transport is affected by interactions of the ionized layer with neutral particles. Nitric oxide (NOx) transport is modulated by plasma behavior (Knipp et al., 2017), and in turn, NOx modulates the warming and cooling rates. Material in the lower thermosphere, when transported, can manifest in the mesosphere as well, for example, in the formation of polar mesospheric clouds (PMCs) (Stevens et al., 2003).

However, winds and transport in the lower thermosphere are challenging to investigate observationally. Larsen (2002) aggregated decades of chemical release experiments to show significant shear in the altitude range of 100–110 km. A unique opportunity to observe material transport in the lower thermosphere arose toward the end of the U.S. space shuttle program during which time a number of orbiting sensors were also operational. Shortly after launch, the shuttle would deposit about 350 tons of water vapor at about 100–115 km altitude, at an almost identical geographic location every time. This water vapor was then tracked by a combination of satellite and ground‐based observations over subsequent days.

Early work on shuttle plumes tied them to evidence of PMC formation (Stevens et al., 2003). Later Stevens et al. (2005) detected the photodissociated water vapor plume via Thermosphere, Ionosphere, Mesosphere, Energetics, and Dynamics (TIMED) satellite global ultraviolet imager (GUVI) Lyman‐alpha observations and related it to PMCs in the Antarctic. Siskind et al. (2003) showed the water to be detectable in TIMED Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) water vapor measurements. Water vapor deposited at about 100 km altitude after the final shuttle launch in 2011 from Kennedy Space Center, Florida, was detected in the lower thermosphere over the next several days. Transport over the Atlantic resulted in observed increases in water vapor in the sub‐Arctic and in PMC brightness in the Arctic (Stevens et al., 2012). More recently, an experiment depositing water vapor in the mesosphere showed the vapor locally cooled and increased the frost point, leading to the rapid formation of PMCs (Collins et al., 2021).

Meier et al. (2011) surveyed launches between 2002 and 2007 whose water vapor plumes were detectable with GUVI and SABER measurements. In seeking to explain the cause for the observed global‐scale observations of water vapor at these altitudes, Meier et al. (2010) discussed that molecular diffusion could transport water vapor, Stevens et al. (2014) attributed the observed distributions to advection and diffusion, and Kelley et al. (2009) argued that two‐dimensional turbulence gave rise to an inverse cascade. Yue et al. (2013) investigated the fast zonal transport of the shuttle plume of mission STS‐121, sustained over several days using TIME‐GCM. With the successful launch of the commercial space transport era for both cargo and passengers, the transport of anthropogenic material as well as naturally occurring material is of renewed importance.

A primary tool for investigating lower thermospheric advection is atmospheric modeling. For example Killeen and Roble (1986) investigated neutral parcel transport at about 120 km and the effects of ion drag due to the E layer ionosphere with a thermospheric general circulation model. A number of surface‐to‐space atmospheric models have been developed that have variants that include high‐altitude charge‐neutral interactions and that can assimilate observational data (H.‐L. Liu et al., 2018; J. McCormack et al., 2017). Most commonly atmospheric models are specified in an Earth‐fixed Eulerian frame as a function of longitude, latitude, pressure, and time. The outputs may be analyzed by mapping instantaneous properties, or by considering zonal mean quantities over time. Polar vortices are of particular interest, but the Eulerian definition of vortex strength can be problematic because it is often frame‐dependent. In some cases, it is desirable to consider properties of a flow that are objective, or observer‐independent (Haller, 2015). The Lagrangian frame advects with the fluid, and enables objectivity.

Prior work investigating the possibility of Lagrangian analysis for explaining material transport used the Horizontal Wind Model 2014 (HWM14) (Drob et al., 2015) to simulate the flow during the last shuttle launch. Wang et al. (2017) hypothesized that there was a Lagrangian coherent structure (LCS) in the thermospheric horizontal flow, that is, a barrier that inhibited the water vapor from reaching an Arctic site. The previous investigation used the empirical HWM14 and considered one water vapor plume deposition event.

In this work, we are interested in exploring further whether LCSs, using simulated flow fields that are more realistic for a specific event than HWM14 because they assimilate data about the current conditions, can help to constrain observed transport. An outstanding question in the space shuttle events is that of the apparently strong meridional transport that has been observed. We aim to find whether there are LCSs derived from the models that could indicate that strong meridional transport should be expected over the course of the initial day after deposition. We focus on two forms of the Whole Atmosphere Community Climate Model with thermosphere‐ionosphere eXtension (WACCMX) that incorporate observational conditions in different ways, and five shuttle launch events distributed over all seasons. Since there is some uncertainty in the capabilities of any model, and in light of the fact that neither model is directly constrained at the altitude investigated in this work, it is useful to see how the model uncertainty impacts the LCS.

We run these models for the space shuttle vapor plume events and output the zonal and meridional winds at an altitude within the vertical range of plume deposition. We use the wind fields to compute the LCSs predicted by the model and compare the model estimates of LCS locations and of tracers to the observations of GUVI Lyman‐alpha. We compare different models' estimates of LCS locations to each other. Finally, because the results between the models differ, we compare the LCSs from the models with zonal wavenumber filtering and at different altitudes. Section 2 reviews the concept of Lagrangian analysis and the whole atmosphere models used. Section 3 describes how we compare the models to each other and to the observations. Section 4 shows the LCSs derived from the models for five space shuttle plume deposition events, with discussion. Section 5 provides closing remarks.

Background

Lagrangian Coherent Structures

A LCS is a manifold (ridge or surface) of maximal separation or convergence in a time‐varying flow. The finite‐time Lyapunov exponent (FTLE) is a commonly used means of identifying the LCS (Shadden et al., 2005). The FTLE map is a geometric approach for defining coherent structures that is “simple” and “practical,” among methods for defining Lagrangian coherence (Allshouse & Peacock, 2015). The FTLE is a scalar field that measures the maximal stretch undergone at a given initial location for a given finite time interval. The computation of the FTLE for a two‐dimensional flow, such as those shown in this work, is summarized here.

Let represent position in a two‐dimensional (2D) fluid domain, as shown in Appendix A. Suppose velocity fields over a domain are supplied over time. The flow map is the transformation of any initial fluid element location at time to a final position at time after a finite time interval : where the integration happens over a particle path rather than a fixed Eulerian grid point. The flow map depends on the initial time and the integration duration. Linearization of about an initial point yields the Jacobian of the flow map at a point , which for a 2D flow quantifies how far apart points that are initially infinitesimally close to end up being:

Interested readers may refer to Figure 1 of Ramirez et al. (2019) for an illustration of the terms in Equation 2. For a fluid element at , the amount of stretching it undergoes can be measured by the FTLE, the largest singular value of :

The FTLE is a scalar field over the domain whose highest values correspond to the points of maximum separation over the time interval. The logarithm and normalization by time interval are conventions in the definition following the Lyapunov exponent formulation, in which linearization of a dynamical system yields exponential growth or decay over time. The LCS is the ridge (or surface in a 3D domain) of the locally maximum FTLE. The LCS represents a locus of the points at time that separates regions of the flow. The LCS is frame‐independent because all observers compute the FTLE field by integrating Equation 1 over the particle path. The LCS structure depends on the initial and final times of the flow map . For a given initial time, increasing has the effect of lengthening LCS ridges. Also, the locus of points comprising the LCS at time themselves advect with the flow. The persistence of an LCS at a given location for a given but different initial times indicates a flow is only slowly time‐varying. LCSs do not directly depend on diffusion rates but may implicitly include the effect when the velocity fields provided result from both advection and diffusion.

Figure 1 illustrates the interpretation of the LCS with an example flow of a time‐varying double‐gyre. Two counter‐rotating vortices alternately expand and contract relative to each other over time with a period of 1 time unit. The flow field at one instant is shown in Figure 1a. The FTLE values over the flow domain are plotted in color in Figure 1b for an integration time time unit. We subdivided the integration time into 50 intervals of duration , computed flow fields at each time, and integrated over one complete period to compute the flow map . The LCS is the yellow ridge that runs up the center and curves left near the top. It separates different regions of the initial flow domain. To illustrate the significance of the LCS, particles are initialized equidistantly at , and their subsequent locations over time plotted at regular intervals until they arrive at final points . Initial particles and on the same side of the ridge end at points both in the left half of the domain. Initial points and lie on either side of the LCS ridge. Point ends up at , in the right half of the domain, unlike . The LCS at time separates from , which indicates that they will undergo different flow trajectories over time. The set of points at which the LCS lies at time also advects with the flow over time (not shown in this still image). For this reason, although particle at some time appears to cross over the LCS ridge (plotted at instant ), it does not actually cross; the points of the LCS at time have moved elsewhere by time as well.

Figure 1

Example two‐dimensional flow of a time‐varying double‐gyre (a) flow field at one epoch (b) FTLE map at and integration time unit of time, e.g., a second, for the flow whose initial velocities are plotted in (a) but vary periodically over 1 s. The colorbar indicates the FTLE value in units of . The LCS at time is the yellow ridge defined by the locally maximum FTLE value. The position of three initially equidistant tracers are plotted over time, ending at final locations .

In this work, the ionosphere‐thermosphere algorithm for Lagrangian coherent structures (ITALCS) is used (Wang et al., 2018) to compute forward‐time LCSs, meaning ridges of maximal separation. Velocity fields are provided as inputs on a global geographic latitude and longitude grid, in m/s in the eastward and northward directions. Since the domain is specified in latitude and longitude, velocities are converted to angular rates in deg/s, and particles are wrapped around as described in the appendices of Wang et al. (2018). Flow fields are bilinearly interpolated in space for integrating particle velocity along a path, and they are interpolated in time between flow field cadence as described in Appendix A.

Atmospheric Models

WACCMX is a whole atmosphere model that incorporates the necessary chemical, dynamical, and physical processes to simulate the troposphere, stratosphere, mesosphere, thermosphere, and ionosphere (H.‐L. Liu et al., 2018). WACCMX extends from the surface to hPa (about 500–700 km depending on solar activity). The horizontal resolution is in latitude and longitude, and the vertical resolution is 0.25 scale height above 0.96 hPa (about 50 km). The finite volume dynamical core in WACCMX includes a fourth‐order divergence damping operator (Lauritzen et al., 2012). This is necessary, and typical of numerical models, to reduce small‐scale numerical noise and instabilities, which can lead to model failure. The resulting winds are influenced by the damping, so in general, would impact the LCSs. However, the purpose of damping is to reduce non‐physical, numerical effects which would lead to LCSs that are not representative of the physical atmosphere's behavior. A detailed description and validation of WACCMX can be found in H.‐L. Liu et al. (2018) and J. Liu et al. (2018).

It is necessary to constrain the model meteorology in order to simulate specific time periods and then compare model output to observations, as we do in this work. This can be done using two different approaches. The first approach, referred to as Specified Dynamics WACCMX (SD‐WACCMX), constrains the model dynamics (i.e., temperature and winds) up to ~45 km based on the NASA MERRA 2 reanalysis (Smith et al., 2017). From 45 to 55 km its dependence on the reanalysis decreases to 0 as the relaxation time increases to infinity, and from 55 km upward the solution is free‐running, including in the analysis region in this work. For the present study, we use the output from existing SD‐WACCMX simulations (https://doi.org/10.26024/5b58-nc53).

In the second approach, the model meteorology is constrained directly by observations using the Data Assimilation Research Testbed (DART) (Anderson et al., 2009) ensemble Kalman filter. WACCMX + DART assimilates conventional meteorological observations in the troposphere‐stratosphere, and also assimilates satellite temperature observations from TIMED/SABER and Aura microwave limb sounder (MLS) up to ~100 km (Pedatella et al., 2018, 2014). This provides an additional constraint above 55 km that is not present in SD‐WACCMX, and can improve representation of the dynamics in the mesosphere (Pedatella et al., 2018).

For the WACCMX + DART model runs, the Aura/MLS temperature observations are assimilated up to hPa (~90 km). The vertical resolution of Aura/MLS varies with height from ~3 km in the upper stratosphere to ~13 km at hPa. The TIMED/SABER observations are assimilated up to hPa (~95 km) with a fixed vertical resolution of 2 km. Observational errors are specified in the data assimilation based on published values (Dawkins et al., 2018; Livesey et al., 2020). The mean bias between the Aura/MLS and TIMED/SABER observations is also taken into account. The bias varies with altitude; near 100 km the Aura/MLS temperatures are ~5 K colder than the TIMED/SABER temperatures (Eckermann et al., 2018). The impact of assimilating middle atmosphere observations from Aura/MLS and TIMED/SABER on the thermosphere dynamics has yet to be fully explored. Pedatella et al. (2018) did, however, demonstrate improved ionosphere specification in WACCMX + DART compared to SD‐WACCMX. This suggests that the thermosphere is also improved in WACCMX + DART.

WACCMX + DART is run for an ensemble of 40 members, generated with small perturbations in wind and temperature at the beginning of each simulation. The WACCMX + DART ensemble mean zonal and meridional winds are used for the LCS analysis.

SD‐WACCMX is more computationally efficient, making long term simulations more feasible. SD‐WACCMX simulations have also been performed for long time periods (1980–2017), enabling comparison for multiple events. WACCMX + DART generally better represents the global‐scale dynamics of the mesosphere and lower thermosphere (MLT) compared to SD‐WACCMX due to the assimilation of observations in the MLT. However, the ensemble approach taken in WACCMX + DART is more computationally expensive, making it more feasible to use only for short case studies. For this reason, WACCMX + DART is used for only a subset of the events investigated in this work, but SD‐WACCMX is used for all events.

Method

In this work, we first consider four space shuttle mission (STS) events for which GUVI observations are available. Table 1 summarizes the events and run parameters used to compute LCSs from SD‐WACCMX winds. The dates selected are: April 8, 2002 (northern hemisphere spring day of year (DOY) 98), October 7, 2002 (fall DOY 280), January 16, 2003 (winter DOY 16), and July 4, 2006 (summer DOY 185). Plume deposition occurs about 220 s after launch through the shuttle main engine cutoff at 512 s (Meier et al., 2011). During the period of deposition of water vapor in the thermosphere the shuttle ground track curves northeastward between ( N, W) and ( N, W) (Stevens et al., 2012). Deposition occurs over a geometric altitude range of about 100 km to about 110 km during a relatively horizontal part of the launch trajectory. Note that STS‐107 had a higher altitude deposition range of up to 115 km, as the only non‐International Space Station (ISS) mission of the flights considered here.

Table 1

Shuttle Plume Events and ITALCS Run Configurations for Comparison to GUVI

Launch	Spring	Fall	Winter	Summer
Year	2002	2002	2003	2006
Date	April 8	October 7	January 16	July 4
DOY	98	280	16	185
STS	110	112	107	121
t_0= Launch UT	20.74	19.76	15.65	18.63
τ= Hours from launch to data a	15.62	18.04	24.44	24.37
DOY of GUVI data	99	281	17	186
Target geometric height (km)	109	109	109	109
Geopotential height (km) for day 1, day 2	108	107	107, 108	108
pressure (hPa)	8.47e‐5	1.09e‐4	8.47e‐5	6.60e‐5

Meier et al. (2011).

The GUVI Lyman‐alpha observations subsequent to each of these launches are mapped. The data is processed following the method in Meier et al. (2011). Superpixels containing 14 raw pixels are created from the along‐track pixels and across‐track pixels of the Lyman‐alpha GUVI “color” (see e.g., Paxton et al., 2004). Next, the Lyman‐alpha anomaly is found by subtracting a scene adjacent to the plume, which removes the high background and viewing‐geometry‐dependent factors. The result is a radiance anomaly in units of R.

We compare maps of GUVI Lyman‐alpha anomaly to the LCS ridges derived from SD‐WACCMX model runs. Based on the fractional depositions as a function of altitude computed by Meier et al. (2011) for STS‐107 (ranging from 100 to 115 km altitude) and for the ISS missions (from 100 to 110 km), balanced by a desire to analyze a similar altitude across events, we select a single geometric height (gmh) of 109 km for analysis of all four events. This is a simplification since the start and endpoints of deposition will be at different heights even within a single launch (see Figures 2 and 3 of Meier et al., 2011) but is consistent with the approaches of other studies (e.g., Yue et al., 2013 used 105 km altitude). As the LCS analysis performed here is 2D, no vertical motion is considered. For each event, the model pressure level that most nearly corresponds to a geometric height of 109 km for grid points between the latitude range of to N and longitude range of to W is used. This pressure level corresponds to a geopotential height within 1–2 km of the 109 km geometric height as shown in Table 1. SD‐WACCMX data are provided tri‐hourly.

The transport calculation for computing the FTLE map includes zonal (U) and meridional (V) wind components but no vertical transport. The water vapor is treated as a passive tracer. The WACCMX models include horizontal diffusion, so the LCSs that result implicitly include diffusive processes that alter the winds. Chemical diffusion that does not impact the winds will be absent from the LCS interpretation. However, previous studies (Yue et al., 2013) showed fast transport due to advection alone, to which the LCS calculation is comparable.

ITALCS outputs the FTLE map and traces individual particles initialized at user‐specified locations through the flow for visualization purposes. Since the shuttle launches from the same position on Earth, the latitude and longitude at which it begins deposition of water vapor in the thermosphere is fairly constant near point = ( N, W). The end point of water vapor deposition is fixed as = ( N, W). We initialize tracers and at positions and , respectively, and track them based on the model flows until they arrive at final positions . Forward‐time LCSs, which indicate regions of flow separation, are computed in this work. Backward time LCSs show regions of convergence but are not plotted in this work for simplicity.

Launch time is chosen as in ITALCS because it is a single well‐defined event only a few minutes prior to water vapor deposition, whose time difference from the start of deposition is much smaller than the model output cadence. The GUVI observations are made at some elapsed time after launch, which is selected to be integration time . The FTLE map is computed by integrating SD‐WACCMX flow fields from the model output from to for particles initialized at each WACCMX grid point. Because the integration time differs for each event and is used in the normalization of the FTLE, the absolute magnitudes of the resulting FTLEs will differ between events. For each event, we are interested in identifying the locus of locally maximum FTLE as the LCS ridges. We also plot tracers of the initial plume deposition; subsequent locations of tracers ending at illustrate the horizontal transport of the water vapor, according to the model.

We then perform inter‐model comparisons between SD‐WACCMX and WACCMX + DART. Table 2 shows the configurations of the different inter‐model comparisons conducted. The first comparison examines how much WACCMX + DART LCSs differ from those of SD‐WACCMX for the July 4, 2006 event (which was also listed in Table 1 because GUVI observations are available) at 109 km. We also compare LCSs as produced by SD‐WACCMX and WACCMX + DART to the GUVI observational data for this event. We select the July 4, 2006 shuttle event because the other events in Table 1 predate the Aura satellite, which launched in July 2004, and whose data improve the quality of WACCMX + DART fields. Note that while the WACCMX + DART run has hourly temporal cadence, the model outputs are downsampled to tri‐hourly to align with the SD‐WACCMX temporal resolution.

Table 2

Shuttle Plume Event and Analysis Configurations for Comparison of SD‐WACCMX to WACCMX + DART

Comparison	1	2	3
Year	2006	2010	2010
Date	July 4	February 8	February 8
DOY	185	39	39
STS	121	130	130
t_0= Launch UT	18.63	9.23	9.23
τ (hours)	24.37	24	24
GUVI Lyman‐alpha data available?	yes	no	n/a
Target geometric height (km)	109	109	60
Geopotential height in SD‐WACCMX (km)	108	108	61
Geopotential height in WACCMX + DART (km)	107	107	59
Pressure in SD‐WACCMX (hPa)	6.60e‐5	6.60e‐5	0.162
Pressure in WACCMX + DART (hPa)	6.60e‐5	6.60e‐5	0.209

The second comparison looks at SD‐WACCMX and WACCMX + DART LCSs for an event on February 8, 2010 at 109 km. The event is chosen based on the list of shuttle missions (Space Shuttle Launches, 2011) for which we have SD‐WACCMX and WACCMX + DART model outputs already available. The GUVI scan mode used for producing Lyman‐alpha anomaly observations stopped working in 2007, so for this event, there are no GUVI observations. However, for this event, the goal is not to evaluate the models using observations but rather to test the range of variability between different model wind fields. Finally, for the February 8, 2010 event, LCSs are compared at a geometric height of 60 km, near the lowest altitude at which SD‐WACCMX becomes free‐running. For the February 8, 2010 event, we also compare the LCSs between the models after filtering the winds as described in Appendix B.

Results

Figure 2a shows the FTLE map for April 8, 2002, at the pressure level nearest to 109 km geometric altitude, with 15.62 h of integration (i.e., the time from launch until GUVI observations are made). The LCSs are the yellow ridges and, at time , mark the points that will undergo the most separation.

Figure 2

Left: SD‐WACCMX FTLE maps at 109 km for events in (a) spring 2002, (c) fall 2002, and (e) winter 2003. Tracers mark the initial plume site and are plotted tri‐hourly as (magenta) and (red) advect, ending at . Right: Maps of only the LCSs at replotted (same color scale as left) and tracers , overlaid on GUVI Lyman‐alpha anomaly column emission rates in R, in grayscale, at for (b) spring, (d) fall, and (f) winter.

Tracer in magenta marks the start of water vapor deposition and is traced over time to the final position , with intermediate points plotted tri‐hourly. Tracer in red marks the site at which the initial water vapor deposition during launch ends, and is also tracked over time to point . Both tracers traverse a circuit clockwise during this time.

Water vapor deposition occurs on the points along the line segment between and . There are LCSs to the south and to the east of the deposition site. Since and are on the same side of each LCS, according to the model they do not separate from each other much over the time interval considered. Figure 2b replots the tracers and the LCSs, extracted as the points whose FTLE exceeds 2.5e‐5 . These are overplotted on GUVI observations of Lyman‐alpha, in grayscale. The tracers appear to end up slightly north of the peak Lyman‐alpha observations. Even though the model's individual tracers do not end up where the GUVI Lyman‐alpha observations appear, the presence of an LCS to one side of the deposition path indicates that the water vapor‐deposited along the ground track is predicted by the model to undergo little stretching. The altitude chosen for the simulation will affect the derived stretching, as the plume deposition occurs over a range of altitudes. The GUVI observations are compact, spanning about in longitude and in latitude. Note that the LCSs are plotted for the locations at , whereas the GUVI data are from scans at time . The fact that the GUVI image of Lyman‐alpha column emission largely falls within the LCSs is, we believe, a coincidence because the LCS ridges mark locations in the flow at time (not ) that separate different flow regions. The locus of LCS points will in general also advect in the flow over time.

Figure 2c shows the FTLE map for October 7, 2002 after 18.04 h of integration. SD‐WACCMX indicates a strong LCS ridge oriented northwest‐southeast seeming to pass through , and a slightly weaker LCS that originates at appearing to weaken extending southwestward. These LCSs from SD‐WACCMX indicate some separation of from , with the magenta tracer traveling southward over the Caribbean and the red tracer traveling north, ending at N latitude over Massachusetts. The tracer points and LCSs are plotted on the GUVI observation map in Figure 2d. In this event, tracer ends beyond the observations of Lyman‐alpha (as with the previous event shown in Figure 2b). Moreover, the LCSs indicate that significant spreading of the water vapor plume should occur. The observation spans about of latitude and of longitude. In the event of Figures 2a and 2b, both are on the same side of an LCS, and the observations of Lyman‐alpha are contained within a narrow latitudinal and longitudinal extent. In Figures 2c and 2d, there are again LCSs but they appear to separate the tracers in the model. In reality, the observations span but the observations are contained to the south and west of the deposition site as though there was not significant spreading of the plume. In this case, the location of the LCS and the northeastward movement of tracer do not agree with the observations. Note that the LCS is the set of points at that undergo maximum separation. This set of points will themselves advect as time elapses, but as they advect, the points do not necessarily remain the points of maximal stretch.

Figure 2e shows the FTLE map for January 16, 2003 after hours of integration. Because the integration time is longer than for the previous event whose hours, the LCSs appear to have a different strength, but this is simply a result of the time‐normalization in Equation 3. Here, SD‐WACCMX modeling indicates a weak LCS ridge (light green) between and . Tracer travels eastward, while travels east and turns south. In this event, plume spreading is predicted by the model and also present in the observations, with up to longitudinal extent (Figure 2f). There is a stronger distant ridge running northeast over northern Canada, but no deposition occurs there.

In summary, we interpret Figure 2 as follows: LCSs immediately to one side of the deposition track indicate containment (spring 2002); LCSs between deposition start and endpoints indicate spreading (winter 2003); if an LCS is predicted in the vicinity of the deposition, it may erroneously predict spreading when containment is observed (fall 2002). Since for the fall 2002 event in Figures 2c and 2d, the water vapor emissions are somewhat contained even though the model's LCS predicts latitudinal spreading, we next examine the sensitivity of LCS location to the choice of model constraint. The October 2002 event was too early for there to be sufficient data with which to constrain WACCMX + DART, so instead, we compare the LCSs for the July 4, 2006 shuttle plume event as predicted by SD‐WACCMX and by WACCMX + DART to GUVI Lyman‐alpha observations.

Figure 3a shows the SD‐WACCMX FTLE map for July 4, 2006 after 24.37 h of integration. There are no significant LCSs in the vicinity of the water vapor plume deposition between and . Weaker LCS ridges are identified west and farther north. The tracers of the start and endpoints of the plume at and both circulate clockwise, with completing almost a complete circuit over the western Atlantic and completing a half‐circuit ending off the coast of northwestern South America. Tracer circulating over the Atlantic ends at well within the area of enhanced Lyman‐alpha emissions in the GUVI observations as shown in Figure 3b, though does not appear to.

Figure 3

Left: FTLE maps for the July 4, 2006 event as computed from (a) SD‐WACCMX and (c) WACCMX + DART. Tracers mark the initial plume site and are plotted tri‐hourly as (magenta) and (red) advect, ending at . Right: Maps of GUVI Lyman‐alpha anomaly column emission rates in R at , in grayscale, with only LCSs at replotted (same color scale as left) and tracers , overlaid as computed from (b) SD‐WACCMX and (d) WACCMX + DART.

Previous tracer simulations for a summer plume event (the July 2011 final shuttle launch) conducted with TIME‐GCM (Yue et al., 2013) and with HWM14 (Wang et al., 2017) showed strong zonal transport over four days and over two days, respectively. For the July 4, 2006 event the SD‐WACCMX model shows only modest zonal transport with individual tracers. The SD‐WACCMX LCSs in Figure 3a are more structured than the LCSs found using HWM14 (compare Figure 3a with Figure 3 of Wang et al., 2017). The peak values using SD‐WACCMX of 3.5e‐5 also exceed the Wang et al. (2017) peak FTLE values (about 1.25e‐5 for an integration time of 2 days), even accounting for the doubling of the time in the normalization of the FTLE. Larger FTLE values are expected here since HWM14 is a climatology, and thus will have weaker, less structured winds. A lack of LCS barriers between does not preclude transport from occurring, but only indicates that there is little spreading predicted in the transport that does occur during that time interval. The area extent of the spread, of longitude and more than latitude, is comparable to or greater than the winter event of Figures 2e and 2f, for which there were at most weak (light green) LCS ridges.

Figure 3c shows the WACCMX + DART FTLE map for the same July 4, 2006 shuttle event after the same duration of integration. In this case, there is an LCS ridge running northwest‐southeast, to one side of the LCS deposition site. The tracers undergo very little separation (the FTLE values are ). One might expect based on this FTLE map that the water vapor would undergo little spreading, yielding vapor plume spread akin to that of Figure 2b, of longitude and latitude. In Figure 3d, the tracers do fall within the GUVI Lyman‐alpha anomaly observations, but the observations extend well north and east of the region spanned by the tracers. The model FTLE map gives no indication that we might expect this much plume spreading.

In summary, the inter‐model comparison for the summer case shows that FTLE maps and tracer trajectories can look different, even for a single model using different forms of data constraint. Note the differences in the FTLE maps of Figures 3a and 3c, in which SD‐WACCMX shows more smaller‐scale structures than WACCMX + DART does, not only at the plume deposition site but across the northern mid‐latitudes. The FTLE map may be more structured in SD‐WACCMX since it only represents a single instance, in contrast to the ensemble‐averaged wind field of WACCMX + DART.

We cannot comment on whether the lack of LCSs northeast of the launch site is typical. Gaining a greater understanding of the frequency of LCS occurrence, or lack thereof, over specific regions is the subject of future work. The lack of strong LCSs over the North Atlantic at the instant shown indicates nothing about whether the water vapor between may be transported to the Arctic at later times. The FTLE map only indicates particle stretching for fluid located at time at the position of the FTLE value shown. Other altitudes in the deposition range 100–115 km could very well present a different FTLE map and, given the presence of vertical wind shears at these altitudes, likely would.

Figure 4 shows a snapshot of the flow fields, spatially downsampled for legibility, for July 4, 2006 at 18:00 UT as output by (a) SD‐WACCMX and (c) WACCMX + DART. Spatial differences in the wind directions between SD‐WACCMX and WACCMX + DART are visible, for example, over the North Atlantic around N, E). Over the eastern U.S., the wind speeds are greater in magnitude in WACCMX + DART than in SD‐WACCMX with opposite zonal component. The flow is similar between the models south of Mexico at ( N, E).

Figure 4

Left half: horizontal wind field for July 4, 2006, 18:00 UT at 109 km, corresponding to “Comparison 1” in Table 2, for (a) SD‐WACCMX (c) WACCMX + DART. Wind profiles in the (e) zonal and (f) meridional directions, with mean (solid) and standard deviation (dashed) for points within of latitude and of longitude of the midpoint ( N, W) of the water vapor deposition site, over the time interval July 4, 2006 18:00 UT to July 5, 2006 18:00 UT. Right half: horizontal wind field for February 8, 2010 09:00 UT at 109 km, corresponding to “Comparison 2” in Table 2, for (b) SD‐WACCMX (d) WACCMX + DART. Wind profiles in the (g) zonal and (h) meridional directions, with mean (solid) and standard deviation (dashed) for points within of latitude and of longitude of ( N, W), over the time interval February 8, 2010 09:00 UT to February 9, 2010 09:00 UT.

Figures 4e and 4f compare the wind profiles for July 4, 2006 as a function of geopotential height. Solid lines indicate the mean speed at a given altitude for grid points within of latitude and of longitude of the midpoint of the deposition track ( N, W), and over the 24 h after launch. Dashed lines indicate the sample standard deviation envelope. Blue lines correspond to SD‐WACCMX and red to WACCMX + DART for (e) zonal winds and (f) meridional winds. The zonal wind mean and sigma envelopes for both models are similar to each other throughout the altitude range. The values vary over space and time, as indicated by the dashed lines. The standard deviation increases with height in both models, indicating greater spatial and temporal variability in the instantaneous winds at higher altitudes. At 107 km geopotential height (corresponding to 109 km gmh), both models show an average eastward zonal wind of 10–20 m/s with a standard deviation of about 30 m/s. The mean meridional winds differ most between 70 and 120 km, but at the geopotential height of 107 km investigated are comparable at about 10–20 m/s southward. The meridional standard deviations also grow with altitude, and at 107 km are about 40 m/s for WACCMX + DART and about 30 m/s for SD‐WACCMX.

The July event had similar means and standard deviations of winds between the two models at the assumed water vapor plume deposition height (107 km geopotential height), and yet the two models yielded somewhat different FTLE maps and tracer trajectories. We compare the modeled LCSs for one final space shuttle event, February 8, 2010. By this time, the GUVI scan mode was no longer operational, so there are no GUVI Lyman‐alpha anomaly maps with which to compare. However, we conduct an inter‐model comparison to analyze the extent of difference between FTLE maps derived from the two models SD‐WACCMX and WACCMX + DART.

Figure 4 plots flow fields for February 8, 2010 at 09:00 UT as output by (b) SD‐WACCMX and (d) WACCMX + DART. Large scale similarities exist between the models, for example, the North Atlantic. Spatial differences in the wind directions between SD‐WACCMX and WACCMX + DART are also visible, for example, around ( N, E). Over the eastern U.S., the wind speeds are greater in magnitude in WACCMX + DART than in SD‐WACCMX, as with the July 4, 2006 event.

Analysis of the simulation results for both the July 4, 2006 and the February 8, 2010 events (not shown) demonstrates that the longitudinal structure of the wind field at these altitudes is dominated by the migrating semidiurnal tide, which has amplitudes of 30–40 m/s in the meridional and zonal wind. The diurnal tide, which has amplitudes of 15–20 m/s at these altitudes, has a smaller contribution to the net wind field. As discussed by Siskind et al. (2003), the tidal phase combined with the local time of the shuttle launch can influence the water vapor transport, making accurate simulations of the tidal variability crucial to correctly capture the transport of the shuttle plume. Yue and Liu (2010) and Yue et al. (2013) furthermore argued that the quasi two‐day wave (QTDW), which peaks in the summertime lower thermosphere at mid‐latitudes (Chang et al., 2011; J. P. McCormack et al., 2014), could play an important role in the fast meridional transport of the shuttle plume to the polar summer.

Zonal and meridional wind profiles as a function of geopotential height are shown for February 8, 2010 in Figures 4g and 4h, respectively, with blue lines for SD‐WACCMX and red for WACCMX + DART. In this case, the zonal wind means at about 107 km geopotential height differs between the models by 40 m/s including a directional difference. The mean meridional wind differs by about 20 m/s. The differences between the means taken over both space and time are significantly larger than for July 4, 2006. We anticipate that, with such different means and standard deviations of 50 m/s at 107 km, the FTLE maps will differ even more than they did for the July 4, 2006 event.

The global LCSs near 109 km after 24 h integration time as derived from SD‐WACCMX and WACCMX + DART are shown in Figures 5a and 5c. In SD‐WACCMX both tracers are surrounded by two LCS ridges, one running east‐west across North America to the north of both points, and one aligned east‐west across the Atlantic to the south of them. From the earlier study of Figure 2, this would imply containment of water vapor‐deposited. For WACCMX + DART, the tracers are separated by an east‐west ridge stretching over the Atlantic. This LCS would predict spreading. In each of the cases, the individual particles have significantly different paths and final endpoints. In SD‐WACCMX, both tracers travel southeast with ending on the northern coast of South America and ending over the Caribbean. In WACCMX + DART tracer travels northeastward, ending at a significantly different place (near western Europe) than it does in SD‐WACCMX. Tracer travels due east to the mid‐Atlantic.

Figure 5

LCSs for February 8, 2010 at (a) SD‐WACCMX 109 km geometric height (b) SD‐WACCMX 60 km (c) WACCMX + DART 109 km (d) WACCMX + DART 60 km (e) SD‐WACCMX vs. WACCMX + DART FTLE difference map at 109 km (f) SD‐WACCMX vs. WACCMX + DART FTLE difference map at 60 km.

Based on the tracers' eastward movement in both models, we might expect significant zonal transport. Based on the presence of LCSs near the deposition site, we might expect significant spreading of the plume if using WACCMX + DART for prediction and compactness if using SD‐WACCMX. Even though both models show LCS ridges across the U.S. and across the Atlantic, the differences in locations appear to correspond to individual tracers having very different trajectories from one model to another. The difference in the 109 km LCS ridge location can be seen in Figure 5e, which plots the difference between the FTLE values of SD‐WACCMX and WACCMX + DART. In this plot, the color scale varies between . Green indicates identical FTLEs at a grid point, yellow indicates that SD‐WACCMX has the larger FTLE, and blue indicates WACCMX + DART does.

Across the U.S., the SD‐WACCMX FTLEs exceed WACCMX + DART over the southeastern U.S. (near the plume deposition points), and appears more structured generally, even though the SD‐WACCMX result uses less data than WACCMX + DART. These effects in the FTLE map likely arise from the fact that SD‐WACCMX represents a single run whereas WACCMX + DART is an average over ensemble members. Overall there are large differences in the FTLE values (and hence LCS ridges) over the equatorial zone as well. SD‐WACCMX shows equatorial LCSs but WACCMX + DART does not. Globally the mean difference in FTLE is 1.21e‐6 and standard deviation in FTLE difference of 1.14e‐5 .

We compare the same models' LCSs at 60 km altitude in the lower mesosphere, shown in Figure 5b for SD‐WACCMX and Figure 5d for WACCMX + DART. Each model at 60 km shows a weak North Atlantic LCS, and to some extent a weak east‐west ridge over southern North America. Here, the tracers all have similar movement, eastward over the North Atlantic. While is at almost the same site in both models, SD‐WACCMX and WACCMX + DART show ending off the west coast of North Africa and over North Africa, respectively. Correspondingly, the LCS structuring is more similar between the models at 60 km compared to 100 km. The FTLE difference plot is shown in Figure 5f. There is more green (an FTLE difference of 0) over the entire region. The mean difference in FTLE globally is and standard deviation is , both about half the magnitude at 109 km.

The better agreement in FTLE values, LCS ridges, and tracers at lower altitudes is possibly due to the fact that the MERRA 2 reanalysis constrains SD‐WACCMX dynamics at and below 55 km, and 60 km is closer to this zone. WACCMX + DART, on the other hand, assimilates observations from the troposphere up to 100 km, thus constraining the dynamics around both altitudes shown here. There is also likely better agreement at 60 km because the flow is dominated more by resolved large scales, in contrast to the lower thermosphere, where unresolved, parameterized, small‐scale waves influence the flow.

To test the effect of higher frequency dynamics, we filter out small‐scale spatial variations in each component of the winds by Fourier series fitting up to zonal wavenumber six at each latitude over all longitudes, as described in Appendix B. The FTLE maps for filtered winds are shown in Figure 6 for (a) SD‐WACCMX at 109 km, (b) SD‐WACCMX at 60 km, (c) WACCMX + DART at 109 km, (d) WACCMX + DART at 60 km. The FTLE difference maps resulting from removing high wavenumber variations from zonal and meridional wind components are shown in Figure 6 for (e) 109 km and (f) 60 km.

Figure 6

LCSs for February 8, 2010 filtered by zonal wavenumber for (a) SD‐WACCMX at 109 km gmh (b) SD‐WACCMX 60 km (c) WACCMX + DART 109 km (d) WACCMX + DART 60 km (e) filtered SD‐WACCMX vs. filtered WACCMX + DART FTLE difference map at 109 km (f) filtered SD‐WACCMX vs. filtered WACCMX + DART FTLE difference map at 60 km.

For SD‐WACCMX, filtering velocities alters the tracer trajectories shown in Figure 6a significantly from those of Figure 5a. WACCMX + DART trajectories in Figure 5c vs. Figure 6c are nearly unchanged. At 60 km, filtering produces no significant change to the tracer trajectories in either model. The differences in LCS location at 109 km are smoother and less structured in Figure 6e than in Figure 5e. However, significant differences between the models remain, with SD‐WACCMX predicting LCSs both north and south of the water deposition site where WACCMX + DART does not. Filtering velocities for wavenumber six or lower variations with longitude has no significant effect on the FTLE differences at 60 km (note similarity of Figures 5, 6). A similar comparison between the models for different altitudes and using filtering was conducted for the July 4, 2006 event as well (not shown), and while the degree of the effects was less pronounced than for February 8, 2010, the findings were similar: tracers and LCSs were more similar at 60 km than at 109 km, LCSs differed at 109 km, and filtering did not improve LCS alignment significantly.

Filtering out of small‐scale variations in the wind does not appreciably improve the inter‐model agreement of LCS locations at 109 km. The lack of difference at 60 km can be attributed to the general lack of small‐scale structures at this altitude. This suggests that smaller‐scale features can introduce variability and uncertainty. Pedatella et al. (2019) discuss the issue of error growth in the MLT. Smith et al. (2017) cover how the lower and middle atmosphere controls the MLT region through the generation of atmospheric gravity waves that grow in amplitude with altitude. These may play a role in the differences between LCSs of different models in the lower thermosphere. Since we have filtered out the small‐scales in Figure 6, it is possible that the remaining differences are largely because SD‐WACCMX is constrained up to 50 km but WACCMX + DART assimilates data up to 100 km and closer to our diagnostic altitude. However, it is not obvious that LCSs derived from WACCMX + DART for the purpose of analyzing water vapor deposition would necessarily be better at transport prediction, given the ambiguity of the July 4, 2006 event. There, the WACCMX + DART tracers were within part of the GUVI observations of anomalous Lyman‐alpha, but the LCS nearby seemed to indicate confinement of that plume akin to the spring event in Figure 2b, which did not take place. It is not definitive as to whether WACCMX + DART performs better than SD‐WACCMX on regional scales for this case since neither model is constrained by observations at this altitude. Given the sparsity of data in the MLT, the performance of WACCMX + DART on regional scales (1,000 km) is less certain, though it will be better at reproducing the larger planetary scales (5,000 km).

An assessment of how the assimilated observations and their errors impact the LCSs in WACCMX + DART or in SD‐WACCMX requires performing controlled experiments in which this can be determined and is beyond the scope of this work. The impact of assimilated observations on LCSs will also depend on proximity of the observations to the specific region of interest for the water vapor plume. We thus cannot state with any certainty how these will impact the LCS structures.

Conclusion

LCSs measure advection objectively by quantifying time‐integrated differences in trajectories over particle paths. LCSs are derived from SD‐WACCMX for a spring, summer, fall, and winter shuttle event for which GUVI observations of Lyman‐alpha are available within about a day from deposition. For the spring 2002 event, LCSs nearby are consistent with GUVI observations. LCSs exist between the water deposition start and endpoints in the fall flow field at mid‐latitudes at about 109 km. In the fall 2002 case, the LCS location appears to indicate spreading of the water vapor plume that is not borne out by observations of Lyman‐alpha column emission rates. In the winter 2003 case, there are no significant LCSs implying separation; the tracers do appear to agree with the spread of water vapor material in the observations. The summer event has no significant LCSs near the deposition site. We find agreement of LCSs of SD‐WACCMX with GUVI in three of the four cases, but comparison of the summer event LCSs between SD‐WACCMX and WACCMX + DART shows differences in the FTLE map as well as in tracer transport.

An inter‐model comparison between SD‐WACCMX and WACCMX + DART for the summer shuttle launch event shows that the FTLE maps differ for the same event, and comparing each with GUVI Lyman‐alpha anomaly observations, it is unclear which model is better at predicting water vapor transport. Neither model is directly constrained by observations at 109 km. For an additional event, SD‐WACCMX and WACCMX + DART LCS maps are compared to examine inter‐model differences. LCS position and structuring differ significantly more at an altitude of 109 km in the lower thermosphere, compared to 60 km in the lower mesosphere. Filtering out small‐scale spatial variations in winds smooths out but does not eliminate the significant LCS differences near the deposition site. LCSs' sensitivity to smaller scales, which are much more difficult to accurately model, may explain the large inter‐model variation in the lower thermosphere. As noted by Smith et al. (2017), model results diverge with increasing altitude.

There is more small‐scale LCS structuring in SD‐WACCMX than for WACCMX + DART. The increased small‐scale structuring in SD‐WACCMX results relative to WACCMX + DART is attributed to the tendency of ensemble averaging in WACCMX + DART to reduce the small‐scale structures since they are less certain and thus occur in different locations in the different ensemble members. The differences in LCSs more generally are a byproduct of background models and methods for constraining the model meteorology. A dearth of observations, specifically temperature and wind observations throughout the MLT, makes it difficult to validate models at regional scales in the lower thermosphere. Ground‐based sources such as lidar and meteor radar are available at specific sites, but as a study of material transport, tracking individual particles of material would be challenging with the sparsely distributed ground sites, particularly as much of the transport in the case of space shuttle plumes was shown to have occurred over the oceans. Ongoing or future satellite missions such as ICON may enable tracking of material transport in the manner of GUVI. Assimilation of ICON or TIMED wind and temperature data at these lower thermospheric altitudes could help constrain the results.

While zonally averaged quantities are useful to quantify differences in latitude and altitude, LCSs aggregate differences in flows that vary in latitude, longitude, and time by examining the effect on material transport. LCS analysis has been applied to investigate the transport of nitric oxide in the polar vortex (Harvey et al., 2021) at 90 km and ionosphere‐thermosphere interactions above 350 km (Wang et al., 2021), and points to a potentially new diagnostic method for analyzing lower thermospheric flow. However, the analysis shown here indicates that the effects of small‐scale waves on material transport may be significant and unconstrained by observations in the lower thermosphere. Thus, LCS analysis for material transport in the lower thermosphere derived from models should be used with caution, as the models themselves may disagree in flow fields significantly enough to yield different resultant transport. It is notable that the 2D LCSs differed significantly even within a given model at different altitudes. This indicates that flow shears with altitude may lead to LCSs at different points even within the shuttle water vapor deposition altitude range. Vertical variations in the wind may play a role in spreading the water vapor. In the future, a full three‐dimensional LCS calculation might show how such ridges are connected into walled surfaces that constrain material.