Intra-Annual Variation of Eddy Diffusion (k zz ) in the MLT, From SABER and SCIAMACHY Atomic Oxygen Climatologies.

Atomic oxygen (O) in the mesosphere and lower thermosphere (MLT) results from a balance between production via photo-dissociation in the lower thermosphere and chemical loss by recombination in the upper mesosphere. The transport of O downward from the lower thermosphere into the mesosphere is preferentially driven by the eddy diffusion process that results from dissipating gravity waves and instabilities. The motivation here is to probe the intra-annual variability of the eddy diffusion coefficient (k zz ) and eddy velocity in the MLT based on the climatology of the region, initially accomplished by Garcia and Solomon (1985, https://doi.org/10.1029/JD090iD02p03850). In the current study, the intra-annual cycle was divided into 26 two-week periods for each of three zones: the northern hemisphere (NH), southern hemisphere (SH), and equatorial (EQ). Both 16 years of SABER (2002-2018) and 10 years of SCIAMACHY (2002-2012) O density measurements, along with NRLMSIS® 2.0 were used for calculation of atomic oxygen eddy diffusion velocities and fluxes. Our prominent findings include a dominant annual oscillation below 87 km in the NH and SH zones, with a factor of 3-4 variation between winter and summer at 83 km, and a dominant semiannual oscillation at all altitudes in the EQ zone. The measured global average k zz at 96 km lacks the intra-annual variability of upper atmosphere density data deduced by Qian et al. (2009, https://doi.org/10.1029/2008JA013643). The very large seasonal (and hemispherical) variations in k zz and O densities are important to separate and isolate in satellite analysis and to incorporate in MLT models.

Introduction

The intra‐annual variation of turbulent processes in the mesosphere and lower thermosphere (MLT) are important to the distribution of constituents both within and above this region, including atomic oxygen. Garcia and Solomon (1985) studied these very processes, with findings that have stood the test of time. Upward coupling of gravity waves from the lower atmosphere plays a key role in the vertical mixing and constituent distribution in the MLT. Turbulence is due to convective and shear instabilities, combined with dissipating gravity waves (Fritts & Alexander, 2003; Hines, 1960; Lübken, 1997), account for the majority of disturbances that result in eddy diffusion of constituents with respect to the background atmosphere (Becker & von Savigny, 2010; Gardner, 2018; Swenson et al., 2018). Eddy diffusion constituent transport effects in the thermosphere‐ionosphere‐electrodynamics general circulation model by introducing a global average, intra‐annual variable eddy diffusion coefficient at 97 km to couple MLT oxygen densities with the thermosphere (Qian et al., 2009, 2013). Another approach has been to propagate tropospheric originating gravity waves from below, via a mechanistic model of MLT composition (Becker & von Savigny, 2010; Becker et al., 2020; Grygalashvyly et al., 2012), and also by whole atmosphere community climate model (Garcia et al., 2007; Liu et al., 2018). Filtering gravity waves from below by stratospheric and mesospheric winds is an important aspect of the coupling. Parameterizing subgrid‐scale phenomena in general circulation models is a difficult task, but finer grid resolution in the models of the middle and upper atmosphere are leading to improved representations of GWs in the MLT and overlying thermosphere (see Liu et al., 2018).

The climatological distributions of constituents in the MLT are influenced by the eddy transport processes which redistribute constituents both horizontally and vertically. One of the simplest to understand is the atomic oxygen distribution, initially characterized by Colegrove et al. (1965). Atomic oxygen is produced above ∼100 km by photo‐dissociation of O2, and diffuses downward by eddy processes including turbulence and dissipating gravity waves (see Figure 1 schematics). The method of parameterized eddy diffusion velocity is determined by the loss chemistry of atomic oxygen, which recombines near 87 km. Swenson et al. (2018, 2019), hereafter S18 and S19, respectively, refined the method of Colegrove et al. (1966) to determine the global mean parameterized coefficient profile k in the MLT, using measurements of OH airglow emissions from thermosphere ionosphere mesosphere energy and dynamics (TIMED) SABER (Mlynczak et al., 2013; Russell et al., 1999) over 16 years, as well as from the Envisat SCIAMACHY (Kaufmann et al., 2014; Zhu & Kaufmann, 2018) measurements of both OH and O(1S) over 10 years.

Figure 1

Mesosphere and lower thermosphere O transport schematic. Atomic oxygen is produced via photo‐dissociation of O2 and diffuses downward by both molecular and eddy processes to the mesosphere, where loss occurs through recombination (e.g., Colegrove et al., 1965, 1966).

This study is Part II of the S19 study, in which we extend our k determination and analysis to examine intra‐annual variations (IAVs) within three latitudinal zones: the northern hemisphere low‐to‐mid latitudes (NH, 15–55°), the southern hemisphere low‐to‐mid latitudes (SH, −15 to −55°), and equatorial latitudes (EQ, ±15°). The investigations are being implemented in a sequence, the intra‐annual variation of bi‐weekly zonal averages (IAVs), the global IAV (by area weighting the three zones), which are deviations from the global mean (see S19), and deserves a dedicated discussion. Intra‐annual variations have been studied by Salinas et al. (2016) using SABER CO2 measurements. Variations associated with thermospheric waves and advection have been described by Jones et al. (2014,  2017, 2018) where thermospheric O densities vary with an AO, we surmise to be driven by k . This study specifically focuses on the IAV of k , and additionally, the IAV of the MLT oxygen density, a parameter that varies separately with respect to the determination of the diffusive flux of O, in the MLT. We feel it is important to establish these basic coupling processes and that they are incorporated into models so that more complex issues of advection and circulation effects can be better analyzed and understood.

Method Summary and Discussion

The primary transport mechanism for O is diffusion, where the total diffusive flux (nv), and the diffusion velocity are the sum of the molecular and eddy components (Equation 1 below, see S19). The integral loss rate of O, via chemistry, is assumed to be supplied by the downward diffusive flux (Equation 2).

The method for the determination of k is as follows. Equations 1 and 2 both describe the downward flux of atomic oxygen, where (a) is traditional composition relationships and (b) is driven by the O loss due to chemistry. The chemical processes are described in our previous two studies, S18 and S19. The chemistry in Equation 2 is described in S18, and the rate coefficients (k1, k4, and k6) are from Sander et al. (2011). We begin the analysis by calculating the downward flux with Equation 2. The flux is then divided by the oxygen density for the determination of the total diffusion velocity versus altitude. Using Equation 3, the eddy diffusion velocity is determined by subtracting the molecular diffusion velocity. Finally, Equation 4, a variable component of Equation 1, is used to determine k . k is the parameterized eddy diffusion coefficient which represents the transport due to mixing from dissipating and breaking waves and instabilities.

The vertical eddy velocity is a function of the total density gradient imposed by the scale height, the atomic oxygen gradient, and the temperature gradient. On average, the O loss chemistry drives the slope of the oxygen density and bottom side (below 96 km) O profile. Vertical advection changes the vertical gradient of O, so k is sensitive to vertical advection. Especially note the oxygen dependence in Equation 4 is via the gradient of O, which is insensitive to a change in the oxygen density, given changes in O density affect the numerator and denominator equally, resulting in a null effect on the k . The previous studies (S19, S19) employing this methodology involved global, long term averaging of O profiles which also negate horizontal and vertical advection effects. In this study, we analyze and describe the uncertainties associated with advection processes on the results presented herein, both in Section 5.2) and in Appendix D. The term definitions for the equations are described in Appendix A.

Equation 2 used to calculate the integral loss rate, from 80 km to z, where z is the altitude for which a diffusion velocity is calculated, was similar to that described by S19, with one exception. S19 defined the OH loss process using: where the global mean values of ozone and hydrogen density were used to determine the O loss rate. In this study, that expression was replaced with the reaction: where M is molecular density, N2 + O2. Equation 6 is the unique source of the ozone in Equation 5, enabling the study to directly incorporate SABER and SCIAMACHY atomic oxygen effects on k . A detailed discussion and references regarding the O recombination (Equation 6) result in ozone formation, and at nighttime results in the excited OH emission was provided by these authors in Section 2.1 of S19. The O density and reaction coefficients for the loss rate are unique for each of the 78 temporal/spatial elements. The second and third terms in Equation 2 are described in Appendix A, where the definition of altitude is geometric.

The integral flux in Equation 2 is an upper limit since some O is produced via pre‐dissociation in the mesosphere by the Schumann‐Runge bands (e.g., Frederick & Hudson, 1980). We computed the average hemispherical production rate from the Schumann‐Runge bands (S‐R bands) in the 85–92 km altitude region in Figure 7 of Koppers and Murtagh (1996) to be 8.5 × 104 cm−3s−1 for the overhead sun, and the average nighttime, the hemispherical O loss rate is 1.5 × 106 cm−3s−1. The ratio is 5.6% for the overhead sun, but the average dayside production rate would be or ∼1/2 this value. We are performing a detailed study of O production and loss continuity that will refine this fraction, but the relative intra‐annual variabilities of k are unaffected. The k values calculated herein, are an upper limit where the values are less by the fraction of O produced locally by the S‐R bands. This fraction is comparable to the fraction determined from the Frederick and Hudson (1980) model values used by S19, where a comparable fraction was calculated.

Data and Analysis

The effects of tidal and planetary waves are minimized by averaging. The sample averages of both SABER and SCIAMACHY data were chosen to be 14 days for 10 years of SCIAMACHY data and 16 years of SABER data. The spatial zones are large which minimizes coupling by advection from adjacent spatial/temporal domains.

The three zones (NH, SH, and EQ) were combined with the 26 two week periods to describe the intra‐annual cycle, a total of 78 temporal/spatial elements for analysis. In addition to the observed SABER and SCIAMACHY O profiles, N2, O2, H, and neutral temperatures must be determined for each of the 78 elements. N2, O2, H, and neutral temperatures were computed using NRLMSIS 2.0, MSIS2.0 (Emmert et al., 2021), for each of the 78 elements. The model has recently been upgraded to include a large amount of available satellite and ground based data over the past 20 years. Both SCIAMACHY (2002–2012) and SABER (2002–2018) covered a solar cycle, whose minimum was 2008. We chose an average F10.7 value for the 2002–2012 as a mean value for the model computations. Note that over the solar cycle, there is very little variation of the constituents or temperature at the MLT altitudes (80–105 km), where the k values are being computed. For equatorial conditions, MSIS determination of changes of 100 F10.7 units resulted in 1% change in N2 density and T, and 11% change in O density. The change in O density is consistent with the change in the FUV flux associated with the Schumann‐Runge spectral region and the O2 pre‐dissociation rate (see Lednyts'kyy et al., 2017). The S18 study found negligible inter‐annual variation in k with the exception of the QBO (Quasi‐biannual oscillation) in the EQ region. The NRLMSIS 2.0 model calculations were made for the respective day of the year for a given period, and the spatial location chosen was −40° latitude for the SH, 0° latitude for the EQ, and +40° latitude for the NH, where the chosen values were representative for the zones.

There is a significant spatial and temporal pattern in the SCIAMACHY data base. The sampling frequency versus latitude and time is plotted in Figure 2. There was some sampling bias in the first 2 years of operation (2002–2004) where sampling at all eligible times was under‐sampled relative to later years. The missing data after 2004 is mainly owing to (a) solar irradiance measurements, (b) spectral calibration, (c) relative radiometric calibration, and (d) dark current measurements. In addition to the above calibration measurements, quote ”The ENVISAT orbit solar occultation was restricted to latitudes between 65°N and 90°N. The lunar occultation was performed from half moon to full moon. For periods of 5–8 days per month lunar occultation measurements provided latitudinal coverage from 30 to 90S. The solar scanning strategy is similar to the SAGE II scanning (Mauldin III et al., 1985): during sunrise, SCIAMACHY scans several times over the full solar disc” Bovensmann et al. (1999). SCIAMACHY data is plotted where the signal is large enough the signal‐noise provides a minimal error to the measurement. Note that SCIAMACHY provides data for both OH (80–96 km) and O(1S, 557.7 nm), (88–105 km). In this study, it was required that both emissions were measured for a given 2‐week period. Additionally, the SCIAMACHY data were not used to compute zonally (i.e., NH, SH, or EQ zones) nor global average results, but it is plotted for respective intra‐annual variations within a zone, for relevant seasonal information.

Figure 2

A plot of sampling statistics for SCIAMACHY nighttime measurements, for latitude versus time (2002–2012). The dotted horizontal lines define the limits of the zonal boundaries within which the zonal data were analyzed. There are significant data gaps in the summer for all zones (northern hemisphere, southern hemisphere [SH], equatorial), and at lower latitudes (<−40°) in the SH.

The TIMED satellite inclination is 74°. The satellite was maneuvered through a yaw cycle every 60 days, at approximately the same day each year, to orient the SABER instrument to view in the anti‐sunward direction. The intent was for SABER to yaw on the same days each year, but over time, the satellite altitude has dropped at the rate of about 1 km per year. The inclination of the orbit has not changed. However, the effect of the altitude decrease is for the yaw dates to creep earlier than their original dates. For example, the first yaw of the year used to occur on January 22. That same yaw now occurs in late December. The yaw maneuver, as well as the TIMED orbit geometry, enabled the SABER instrument to acquire data on all days of the year, unlike the SCIAMACHY data described in the previous paragraph. The latitudinal coverage has a sampling bias, sampling further southward in a given cycle, and alternately northward bias in the subsequent cycle. The number of measurements at all latitudes between ±55° is large for all years for the 2 week sampling performed herein. Within a given 60 days yaw orientation, there is a variation in the local time of night sampled through the cycle. The amplitude of O density variation is discussed in the data analysis section that follows.

Equation 1 lacks advection terms (Gardner, 2018), an assumption implying turbulence and wave processes uniquely and solely describe the vertical O distribution. The vertical distribution of O is driven by the downward flux of O, via the diffusion velocity (k ), as well as advection. Our prior analysis (S18, S19) involved calculations of globally averaged k . In those studies, the global average constituent profiles for all latitudes were averaged for a minimum of a year (S18). As a result, any advection contribution was arguably minimized through long‐term averaging. Advection potentially influences the O density from an adjacent zone. It was pointed out in the previous section that k and the O density are treated as separate variables. k (driven by the vertical gradient of O) and O density (where horizontal distributions are potentially influenced by advection) are separate. The separation of variables in the intra‐annual cycle (k and [O]) is critical to establishing the vertical coupling of constituents in the MLT.

The 14‐day interval was chosen with the rationale that follows. Diffusion times are an important consideration. Considering a breaking wave condition in the layer and an anomaly in the altitude distribution is redistributed over altitude by diffusion in time. Lednyts'kyy et al. (2017) measured the time delay from the solar variation in the 27‐day rotation (and associated photo‐dissociation of O2) to the time the variation appears in the O(1S) emission near 95 km, to be ∼13 days. It is estimated that it takes an additional week to diffusively transport O from 95 to 87 km, the altitude of maximum loss via recombination (S19). Consequently, sampling average composition distributions at a temporal resolution shorter than 2 weeks would potentially, fail to reach an equilibrium condition. This criterion is overstated considering averaging for 16 years. A statistically significant number of measurements is also necessary, clearly evident in Figure 2 for SCIAMACHY. Twenty six 2‐week periods constitute the data elements for the analysis of intra‐annual variability, for each of the three latitudinal zones.

The O density from both the SABER and SCIAMACHY data archives were computed for each of the elements for each year available (16 years for SABER and 10 years for SCIAMACHY) and averaged for all the years. As a result, the k determined for each of the elements represents the climatological mean for that element.

Results

Results for the average time evolution of k and a 2‐D (day of year vs. altitude) variation of amplitude for 16 years of SABER data are illustrated in Figures 3a–3e and Figures 3b–3f, respectively. Tables of values of k for Figures 3b, 3d, and 3f are provided in Appendix C. The SCIAMACHY data were not included in this initial analysis due to the sampling biases described above. The 2‐D plot illustrates the dominant periods being an AO, in both the NH and SH below 87 km, and the SAO at EQ at all altitudes consistent with the latitudinal variability in IAVs of received solar radiation and surface temperature (see Picone et al., 2019), as well as observed and model IAVs in k and middle atmospheric winds by Garcia and Solomon (1985) and Garcia et al. (1997). The amplitude of the EQ SAO in spring is larger than in the fall. The lack of variability with altitude in k for the EQ versus the NH/SH is unexpected. The details of the intra‐annual AO and SAO variations for the respective zones are described in detail in the following subsections.

Figure 3

(a), (c), (e) Variation of k deduced from SABER measurements, with the day of year in the northern hemisphere (NH), equatorial (EQ), and southern hemisphere (SH) latitude bands. Each colored line represents a different altitude in the range of 80–100 km. (b), (d), (f) Variation of k with day of year and altitude is shown for the NH, EQ, and SH latitude bands.

Northern and Southern Hemisphere and the AO

The results for the NH low‐to‐mid‐latitude region are shown for SABER in Figure 4a. The most prominent IAV depicted at the four altitudes (96, 91, 87, and 83 km) is the AO, below 87 km, with a maximum in the middle of the summer, and a minimum in winter. We also note that there is minimal variation in k above 87 km, including the upper limit of SABER data, at 96 km. Figures 4b and 4c are plots of k for both SABER and SCIAMACHY data for 83 and 96 km, respectively. Some SCIAMACHY data is missing in mid summer, but with the data available, the suggestion is the amplitude at 83 km of the summer AO maximum may be slightly less than that for SABER (Figure 4b). Figure 4c illustrates there is little IAV at 96 km in both SABER and SCIAMACHY. There appears to be a minor peak near the summer solstice in both SABER and SCIAMACHY.

Figure 4

k versus DoY for the NH. (a) Values determined from SABER O, for 96, 91, 89, and 83 km. The AO dominates at 83 km, driven by the meteorology from below, with a maximum in the summer. (b) Same as (a) except 83 km, for SABER (solid) and SCIAMACHY (dots). (c) Same as (b) except for 96 km.

The AO of k is the dominant oscillation in both SABER and data for the SH (Figures 5a and 5b) and SCIAMACHY for the winter, in b. The altitude of most variability is below 89 km increasing to the lower limit near 80 km (Figures 3b and 3f). The variabilities are a maximum near the summer solstice and a minimum in the winter (near the solstice). The amplitude of the winter to summer k at 83 km changes by a factor of ∼4 (see 3b). The SCIAMACHY and SABER k are nearly equal for the winter, where SCIAMACHY data was available, with a hint of winter to summer transition at day 110. A clear 180 days phase shift in the AO is evident (Figure 5b), clearly associated with the season.

Figure 5

(a) Same as Figure 4a, except for the SH zone. (b) SABER southern hemisphere (SH) and northern hemisphere (NH), and SCIAMACHY SH and NH at 83 km.

k as a function of altitude calculated from SABER and SCIAMACHY O measurements, for the NH summer maximum versus winter minimum is shown in Figure 6. SCIAMACHY on Envisat is in a 10 a.m./p.m. polar orbit. The sampling pattern is the result of the requirement that the full line‐of‐sight of SCIAMACHY shall be in complete darkness and that some calibration measurements are performed in the southern hemisphere during nighttime, resulting in fewer data points in the SH compared to the NH. The data chosen for this figure was for periods 26 (end of year) and 1 (beginning of year) for the winter profile, and period 12 (early June) for the summer values. SABER winter was the same as SCIAMACHY, but the summer was period 13 and 14 (late June and early July), chosen for its availability at the summer solstice. The k plots for both SABER and SCIAMACHY illustrate similar variations, below 90 km, with maximum to the lower limit of the data at 80 km. The NH is plotted since SCIAMACHY has data for both winter and summer for this comparison of k . The large k at the 80 km limit of sensing from SABER CO2 in S19, was a result following the original analysis by Salinas et al. (2016). The integrated loss of O was integrated for the summer versus winter for a difference of 20% in the flux at 96 km, for those two extremes. The discussion relevant to these changes follows in the next two paragraphs.

Figure 6

k versus altitude for the northern hemisphere summer maximum (at summer solstice) and winter minimum (at winter solstice), for SABER and SCIAMACHY. The variability begins near 90 km, increasing with decreasing altitude. The variation between summer and winter is factor of ∼4X for both O climatologies at 83 km.

Figure 7a is a plot of the intra‐annual variation of k at 96 km for the NH and similarly, Figure 7b for O density. The major variation in the O density is a broad peak near the summer solstice. In Figures 7a and 7b, the dates the satellite performed a yaw maneuver every 60 days, directs the SABER viewing direction to be anti‐sunward with respect to the orbital plane. This yaw oscillation (YO) performed nearly the same day each year, is directly correlated to the same periodic brightness variation in the O density (Figure 7b). In each yaw cycle, the local hour sampled changes from the beginning to the end of the cycle, and consequently, brightness variations associated with local time variation contribute to the cycle in O density. There is variability in k (Figure 7a) also, but not directly correlated since the O density has a negligible effect on k . There is one event marked P2 in O density near DoY 170, and a spike in k correlates with the sharp trough, following P2 in O density, which will be described in the discussion section. This event is not just a YO, but rather an event observed by both SABER and SCIAMACHY.

Figure 7

The intra‐annual variation of (a) k and (b) O density for SABER at 96 km are shown with the vertical double arrow illustrating a fractional amplitude indication. The circles are the days of the satellite yaw events, each year. The relatively large ’spike’ at day 170 in k in (a) and depletion in O density (b), is referred to as the P2 event.

The next step in our analysis is to better understand the continuity and downward O flux in context with the AO and summer enhancement in O at 96 km, and the k enhancement below 87 km. The variability of atomic oxygen flux at 96 km has been calculated for the NH and is plotted in Figure 8. The eddy diffusion velocity was determined for the same method used to determine k in Figure 7a was multiplied by the O density (Figure 7b) for the calculation of the flux for Figure 8. An amplitude arrow of 10% (or minimum to maximum of 20%) is indicated in the figure. There is general consistency with the integral loss of O below 90 km, where there is a maximum at the summer solstice and a significantly lower flux in winter.

Figure 8

The intra‐annual variation of the SABER eddy flux in the northern hemisphere at 96 km. The flux has a maximum in summer, with lower values in fall, and significantly lower in spring. The vertical arrow illustrates the required change between summer and winter solstice necessary to account for the change in the loss of O (and increase in k ) below 87 km, shown for SABER in Figure 6.

Equatorial Region and the SAO

A semiannual oscillation in k is clearly evident at the EQ, with a much smaller AO than at low‐to‐mid latitudes (Figure 9). Note the SAO dominates at all altitudes, with a minimum SAO amplitude at 83 km. The amplitude of k at 87 km and above varies between a summer solstice minimum and spring equinox by a factor of ∼2, and from the summer solstice to the fall equinox by a factor of ∼1.5. An observation in the phase shift in the SAO was noted in the spring when the amplitude is largest. The phase propagates upward near spring equinox from day 65 at 83 km to day 100 at 96 km.

Figure 9

(a) k versus DoY for EQ 96, 91, 87, and 83 km for SABER. The spring equinox is the dominant oscillation at all altitudes. (b) Same as (a) except at 96 km for SABER and SCIAMACHY, and (c), same as (b) except for 83 km.

Zonal k and O Density IAVs

Figure 10 is a plot of the average k versus altitude for the NH, SH, and EQ zones using SABER data. Note the NH and SH profiles are almost identical with altitude. Also noted is the near‐constant distribution with altitude at the EQ zone.

Figure 10

Annual average k versus altitude values for SABER, for the northern hemisphere (NH), southern hemisphere (SH), and equatorial (EQ) zones. The NH and SH profiles are nearly identical. Note the EQ profile is almost constant with altitude. The horizontal bars indicate the standard deviation at 92, 87, and 82 km for the EQ and at 95, 90, and 85 for the NH for the 16 years of SABER measurements.

We hypothesize this is likely due to a difference in the sources contributing to k , both damped gravity waves and instabilities. Details of the rationale are described in the Section 5.

Figure 11 is a series of plots of the IAV associated with SABER O density. Figures 11a–11c are the IAVs of the percentage of O density change with altitude with respect to the global average density profile. Figures 11g and 11h describe the density IAV of O density at 96 and 85 km, respectively. The phase shift of the AO at the NH and SH with the season as well as the large SAO at the EQ region at 85 km are dominant features. It is noted in particular that the amplitude of the O density below 87 km is 180° out of phase with k , suggesting the large values of downward diffusion velocity in the summer depletes the O. On the contrary, in the EQ zone, the fact that the enhanced k (and diffusion velocity) is larger at all altitudes, the larger O densities near the altitude of maximum density (96 km), supplies the O density from above, overcoming the O‐losses at lower altitude. This is a major difference between the influence of k in the mid‐latitudes and the AO effect in both hemispheres, to that of the SAO in the equatorial region.

Figure 11

The intra‐annual variations (IAV) of the percentage change in O density versus altitude, relative to the zonal mean altitude profile for the northern hemisphere (NH) (a), equatorial (EQ) (b), and southern hemisphere (SH) (c). (d) The IAV in SABER O density at 96 km and (e) for 85 km. Note the AO amplitude phase switches from a summer maximum at 96 km, to a winter maximum at 85 km, where the minimal O density coincides with the maximum in k at low altitude (see Figures 3, 4, and 5b). Also note the significant EQ spring equinox amplitude at 85 km, in phase with the magnitude and vertical extend of k (Figures 9a–9c), rather than out‐of‐phase. Note event P2 with a maximum near the upper altitude of the sampling altitudes at 96 km near DoY 170 in the NH (a) and ∼180 days later in the SH, (c).

Late Spring Event (NH, P2)

A large spike in O density and k in the NH occurs in the spring in the SABER data shown in Figures 7b, 11a, and 11c, which is also present in the SCIAMACHY O density. Examination of a large amplitude at 96 km is noted with a peak at DoY 150 and a minimum (trough) at DoY 170 in the O density. Note the large peak in k (Figure 7a) coincides with the trough in the O density (Figure 11b) where it has been highlighted with a dashed line (phase progressing upward). The feature is also apparent in the SH, shifted by 180 days.

The NH O density at 85 km versus DoY is shown in Figure 12. Considering the local time spread of the measurements over the SABER night which for ’night’ consists of local times where the SZA  (solar zenith angle, solid line). It was noted, that when the local time hour intervals are made smaller, the amplitude of the Spring event changes, illustrated in this case, for SZA  (dotted line). The SABER data has a local time bias associated with the yaw periods (60 days), where the local time sampled shifts in local time from the beginning to the end of the yaw period. Studies of tidal influence on the mid‐latitudes have recently been demonstrated by Tian et al. (2021), where meteor wind climatology observed the diurnal variability as well as IAV tidal effects on momentum fluxes associated with damped AGWs. We believe the change in amplitude of this spectral feature is due to filtering by the tides. The SCIAMACHY O data is also plotted as squares, which describes the same Spring event feature, where dashes highlight the specific feature amplitude. The SCIAMACHY data has gaps (shown as straight lines). Note that SCIAMACHY is always sampling at 10 p.m. local time. It is also important to note that the O densities are relatively small at 85 km, whereas the fractional changes due to the AO are large.

Figure 12

Northern hemisphere O density versus DoY, at 85 km. The SABER nighttime data (solid line), SABER SZA (dashed), and SCIAMACHY data (squares) are shown. Both SABER and SCIAMACHY clearly identify the P2 brightness enhancement. The amplitude effects of the P2 event include local time sampling biases, which are discussed in the text. Note the SCIAMACHY data gaps and SABER Yaw dates are indicated.

We hypothesize this Spring event may be a brief, vertical extension of the AO in k . The specific source of the effect at higher altitude could be a result of an increase in the source amplitude in the troposphere and/or a change in the integral dynamical filtering effects which primarily occur in the stratosphere. We also note the slope with an altitude of the O density enhancement, is similar to the slope of the EQ k at the spring equinox (SAO). It is plausible the event is a brief, very large excursion of a single or few years. Given the fact the oscillation appears in both SABER and SCIAMACHY data (Figure 12), those years would necessarily be between 2002 and 2012 when both SABER and SCIAMACHY were operating. We examined the O density for extreme outliers as well as the standard deviation of SABER O density. The O density (1011 cm3 for the NH period 12 prior to 2012 (1.11, 1.64, 1.65) versus the entire sample (1.29, 1.62, 1.53) at 92, 94, and 96 km, respectively. There was no appreciable difference noted for a given year. We conclude that oscillation is an annual geophysical feature.

Global Average k IAV

Global average (i.e., between ±55°) of SABER O density and k as a function of day of the year at 96 km is shown in Figure 13a. The 96 km altitude was chosen because (a) it is the highest altitude of O resulting from OH airglow inversions used by SABER with relevant to global means, (b) is near 97 km chosen by (Qian et al., 2009) and (Salinas et al., 2016) intra‐annual variation of k , studies that are compared to later herein, and (c) it is representative of the altitude of maximum O density contribution to the downward flux. The global average for both O density and k is computed by area weighting each zone by its an effective fractional area (0.364 for NH and SH each, and 0.272 for EQ). The Global average IAV (Figure 7b) is dominated by an SAO, reflecting the large EQ SAO contribution to the average. The intra‐annual values of SABER k are also plotted in Figure 13b where it is compared with that of (Qian et al., 2009).

Figure 13

The global average of O density (a) versus day of year for SABER at 96 km (Solid circles). The northern hemisphere, southern hemisphere, and equinox are also plotted for perspective. There is a negligible difference between winter and summer, although the SAO is clearly dominant. (b) The global average k versus day of year for SABER at 96 km (dotted) and Qian et al. (2009) (solid). The mean k for the two plots is nearly identical, but there are significant differences in the intra‐annual variations.

The global average values of k versus altitude for SABER are plotted in Figure 14 and compared with the values from the study of global means by S19.

Figure 14

k versus altitude for the global average (GA) SABER O from this study compared with the SABER O and CO2 and SCIAMACHY O derived results in S19.

Discussion

The phase and a peak altitude of the AO are consistent with dissipating and breaking GWs, which propagate upward from the lower atmosphere during the eastward phase of the stratospheric circulation. Westward propagating, high‐frequency waves are unfiltered during this eastward phase and propagate freely. This hypothesis is also consistent with the extended increase which begins and ends near spring and fall equinox; the times at which the stratospheric winds reverse. These results are consistent with the analysis by Garcia & Solomon (1985). Their analysis of O3 observations and O shape profiles were key elements of their discoveries, which here are confirmed and refined upon with SABER and SCIAMACHY data. These results are generally consistent with predictions from theory (Hines, 1960). A. Z. Liu (2009) analyzed the annual variation of k from lidar observations at 35°N (Starfire Optical Range), where the IAV exhibits a similar peak in amplitude in summer, but at slightly higher altitudes (∼90 km). Enhanced GW activity at mid‐latitudes was also observed by Gardner et al. (2011) and Gardner (2018). Meteor wind observations at mid‐latitude (∼40°N by Tian et al. (2021), clearly illustrate IAV of zonal momentum fluxes to have variability in altitude and season similar to the k variability described in Figure 3b. This study also demonstrates the significant effect the tidal phases have on diurnal variability.

We computed the difference in the integral loss of O between 80 and 96 km in winter versus summer due to the AO, using the SABER profiles shown in Figure 6. The calculated difference in O loss between summer and winter solstice due to k corresponds to a change of 20% in the downward flux of O at 96 km. That difference should reflect the change in either the diffusion velocity (k ) or O density or a combination of both. The change was 20%, from the winter minimum to the summer maximum, or an oscillation amplitude of 10%. We note there is no change in the global average k at the solstice (Figure 6b). According to the observational evidence from SABER and SCIAMACHY, the IAV in the diffusion velocity at 96 km is minimal throughout the annual cycle. The evidence lies in the variation in the O density, and an initial study by Smith et al. (2010). A study of O variation with season was accomplished by Chen et al. (2019), who measured O density oscillations using the GOMOS instrument, and analyzed the AO, SAO, and QBO amplitudes for a few years of observations. These observations complemented studies by Zhu et al. (2015), Lednyts'kyy et al. (2017), and followed by Chen et al. (2019) where amplitudes of the AO were 11%, 7%, and 9.6%, and for the SAO 15%, 12%, and 18%, respectively. The study performed by Chen et al. (2019) had three zones: 20–30°N, −20 to −30°S and an equatorial band. The intra‐annual variation of the O density for SABER NH, SH, and EQ versus day number is shown in Figure 7a, and for both SABER and SCIAMACHY in the NH in Figure 7b. Clearly, the amplitudes of the AO for atomic oxygen are consistent with the climatologically determined loss of O in the MLT.

The IAVs in k at equatorial latitudes exhibit a more prominent SAO (Figure 9), with larger amplitude at the spring equinox (2X) with respect to a summer minimum than the fall (1.5X). The SAO amplitude is reduced below 85 km. The EQ region ±15° is dominated by the influence of the Inter‐Tropical Convergence Zone (ITCZ) in the lower atmosphere, a key factor in forcing the diurnal tides. These results are consistent with the theory described by Dunkerton (1982), which hypothesized that an observed SAO variability in the zonal wind at equatorial latitudes (Hirota, 1978) combined with Kelvin waves selectively enabled gravity waves to propagate into the mesosphere. The enhanced values of k that extend well into the upper mesosphere are consistent with this hypothesis. IAV observations of meteor radar winds at Jicamarca (∼12° N, Guo & Lehmacher [2009]), illustrate strong tidal oscillations, with the largest amplitudes at spring and fall equinox, nearly identical to the equatorial IAV amplitudes of k derived herein. This result strongly supports the consideration that wave‐tide coupling is directly responsible for the larger k amplitudes at the EQ versus NH and SH, in the 80–90 km altitude region (Li et al., 2005). illustrates a form of wave‐tide coupling interaction with the diurnal tide observations at low latitudes, where Mesospheric Inversion Layers associated with vertical mixing and turbulence, form with a tidal phase. In addition, a secondary consideration involves the Quasi‐Biannual Oscillation (QBO). Swenson et al. (2018) described a QBO variation in k at EQ latitudes, also with reduced amplitude at 83 km compared to higher altitudes, likely due to the wave filtering by the QBO at lower altitudes.

The AO in k affects altitudes below 87 km (Figures 3b, 3f, and 5b), whereas the SAO extends to 96 km, the upper limit of SABER (Figure 3d and 9a–9c). The near constant k versus altitude for the SAO enhanced equinox regions is clearly shown in Figures 4 and 5. A possible explanation for the extended altitude region of the SAO is wave‐tide coupling. The large amplitude in the diurnal tide (DT) at equatorial latitudes results in a wave‐tide interaction (e.g., Li et al., 2005). Figure 10a depicts the k profiles for both the annual average for the NH and SH, representing the mid‐latitude and the EQ. The fact that the EQ profile is uniform with altitude, not just at equinox periods but throughout the year, supports the hypothesis that the wave‐tide coupling influences the vertical distribution of turbulence in the EQ region. The fact that the average k is 2–3X larger for the EQ than for the NH/SH is consistent with the ITCZ being a strong convective source of gravity waves. The upward propagating waves experience minimal stratospheric filtering at the equator. Wave coupling with the large amplitude diurnal tide results in a significantly larger eddy diffusion effect in the EQ MLT.

Historically, k (and the diffusion velocity) have been used as a parameter to drive composition effects in a number of general circulation models. Colegrove et al. (1965, 1966) used this approach to define the k relationship to the bi‐directional flux (nv) of atomic and molecular oxygen. S18 modified the original approach and solved for the diffusion velocity only. We understand today that atomic oxygen can be influenced by a host of other considerations including production, loss, and transport by waves, on the scales discussed herein. In a given hemisphere, the summer produces more O via photodissociation than in the winter, whereas the meteorology and forces from below that are responsible for the eddy velocity and downward transport to O loss do not necessarily map to the production timeline. The intra‐annual diffusive coupling between 140 and 96 km is primarily due to molecular diffusion throughout the altitude region, as well as eddy diffusion below ∼105 km. Diffusion upward into the thermosphere reflects the AO that is well documented in the very early thermospheric composition models, for example, Jacchia (1964) where thermal expansion in the thermosphere is important, where O densities vary diurnally by 100%. On the contrary, a 20% change in atomic oxygen density which is the maximum change between the AO solstice extremes, imparts a relatively minimal change to the thermosphere. The k deduced by (Qian et al., 2009), which was derived from thermospheric observations, includes SAO and AO amplitudes where the AO (January and July) amplitude is as large as the SAO amplitude. The total time for atomic oxygen production to the altitude of O(1S) emission was observed by Lednyts'kyy et al. (2017) using the correlation between the emission and solar rotation. The chemical loss of O primarily occurs below 93 km, with a peak loss near 87 km. The O density above 93 km is a reservoir of O that is diffusively and dynamically coupled, wherein there is time variation in that coupling that is dependent on the values of the k and the spectrum of upward propagating large‐scale waves of lower atmospheric origin that dissipate in that region. It is the eddy diffusive process and thus k in models, largely below 93 km, that supports the chemistry of O loss.

Derived k at 96 km is compared to Qian et al. (2009) in Figure 13. The global average k from this study of MLT composition effects retains the EQ zone dominant SAO. The ”top‐down” approach of (Qian et al., 2009) yields a much larger k amplitude, that is out of phase with the dynamical‐chemical balance approach employed herein. The message here is the climatology of the O density and k dominate the vertical transport, and the global mean has little value to a model that is describing the composition since seasonal effects dominate. In particular, the O density and k with the NH versus the SH are dramatically different for a given time, and when analyzing satellite data for a given orbit, the respective hemispherical responses should be kept separate, in order to account for the dominant seasonal influences.

Implications for the Overlying Thermosphere and Ionosphere IAVs

The IAVs in k that we have deduced will naturally have implications for upper thermospheric and ionospheric IAVs (e.g., Qian et al., 2009, 2013; Yue et al., 2019) through seasonally‐dependent transport of atomic oxygen in the MLT region (Jones et al., 2017, 2018). The k analysis by Qian et al. (2009) include both the AO and SAO components, with the AO amplitude at a solstice similar to the SAO amplitude. Figure 13 clearly shows an SAO in global average k that is consistent with CO2‐derived k results from (Salinas et al., 2016), and is much weaker than what Qian et al. (2009) inferred from satellite drag data in the upper thermosphere. Taken with the Salinas et al. (2016), our weaker IAVs in k indicate that the overlying thermosphere and ionosphere SAO is not primarily driven by IAVs in k , but rather acts in concert with the predominant thermospheric spoon mechanism (Fuller‐Rowell, 1998; Jones et al., 2018).

Further, our k deduced from SABER O, is of opposite phase relative to those either produced by gravity wave drag parameterizations or invoked in the NCAR thermospheric general circulation models (see Jones et al., 2017; Qian et al., 2009). This oppositely‐phased IAV in k calculated herein is likely due to IAVs in the SABER O density (see Smith et al., 2010), and thus deduced downward O fluxes calculated for in the MLT region. Differences between IAVs in k deduced from SABER O using Equations (1), (2), (3), (4) and, for example, the NCAR thermosphere‐ionosphere‐mesosphere‐electrodynamics general circulation model is probably because the atomic oxygen flux in Equations 1 and 2 are the total vertical flux of atomic oxygen, including eddy and molecular diffusion and the “bulk” vertical wind (see Jones et al., 2018, 2021). While upper thermospheric general circulation models are able to separate all these different processes, assumptions made in Equations (1), (2), (3), (4) lead to a slight convolution between the eddy diffusive flux and “bulk” vertical wind flux of O.

Potentially, most important for upper thermospheric and ionospheric IAVs are the results presented in Figures 4 and 5, which show a strong AO in k at middle northern and southern latitudes. To our knowledge, unless one uses a large AO in k (like Qian et al., 2009) in the MLT region, upper atmospheric general circulation models do not accurately reproduce the observed AO in thermospheric mass density and ionosphere electron density. Therefore, the latitudinal dependence of IAVs in k within such models should be re‐evaluated given our results. Perhaps, one might expect these upper atmospheric general circulation models would produce a more realistic thermospheric and ionospheric AO in pertinent model parameters if they properly accounted for the latitude and seasonal dependence of k IAVs deduced from SABER. Such reasoning is further supported by recent results from Malhotra et al. (2020), which quantified the sensitivity of middle‐upper thermospheric dynamics, energetics, and composition to changes in O density between the 95–100 km in the global ionosphere thermosphere model (model lower boundary at 95 km).

The temperature profile between 95 km and the thermosphere is an important attribute in the diffusive coupling of composition between atmospheric regions, especially for projecting compositional IAVs into the upper thermosphere and ionosphere (see Jacchia, 1970 and Equations 1 and 4 herein). Temperature measurements from SABER CO2, as well as Gobal‐scale Observations of the Limb and Disk and Ionospheric CONnection Explorer satellite measurements, will provide improved temperature IAVs, for the models above to validate against. The modeling studies above clearly demonstrate that eddy diffusion, neutral wind transport, and temperature all play an important role in MLT coupling of composition with the thermosphere. Further, seasonal production and loss also play a role for O. With these new space‐based assets, providing long‐term datasets, we are starting to be able to truly assess the ability of our general circulation models in the middle and upper atmosphere.”

Effects and Uncertainties Due to Advection

Contamination of the vertical gradient and density of O in a zone can occur as a result of meridional advection from one zone to another, and by vertical advection within a zone. The midlatitude zones (NH, SH), between 80 and 90 km are the zones where there is a strong seasonal difference in k (and diffusion velocity) that are subject to contamination of the O densities by either (or both) vertical and horizontal advection. The AO, driven by the residual circulation below 90 km, is a condition of significant duration (weeks), that has a significant amplitude in both meridional and vertical winds. Analysis for determining these effects on the eddy diffusion coefficients calculated herein involved a close examination of the altitude distributions of the diffusion velocities and time to reach an equilibrium state, the chemical lifetime relevant to advected O from one zone to another, and the climatology of meridional advection velocities associated with contamination from adjacent zones. A detailed analysis is provided in Appendix D supporting the summary given below.

Our findings for the mid‐latitude, AO amplitude of the vertical upwelling at summer solstice was estimated to be 0.4 cms−1 at 80 km. The chemical lifetime at this altitude is ∼1 hr, resulting in penetration of 15 m at that altitude (see Appendix D3), a negligible uncertainty to the vertical distribution of O density and k calculations herein. Similarly, horizontal advection at 85 km, was calculated given a 3 ms−1 horizontal advection velocity at the boundary between the NH and the EQ, with a zone penetration of 10 hr (the chemical lifetime, see Appendix D4). The rationale is that advected O traveling further into the zone would recombine, and its vertical distribution due to diffusion would dominate over advection. The calculated penetration is ∼1° in latitude, contaminating 1/30th (or 3.3%) of the EQ zone. The penetration volume is assumed to be 50% contaminated with respect to chemical loss of O, for uncertainty due to horizontal advection of 1.7%. These relatively small uncertainties allow the latitude zones to be made smaller below 90 km, in future studies.

Above 88–90 km, the diffusion times are less than the chemical recombination time. At 96 km, for example, a 4‐day penetration of 3 ms−1 results in a 10° invasion of the EQ zone resulting in a 33% area contamination, and at 50% contamination, results in a 16.5% uncertainty. In the lower thermosphere the effects and uncertainty due to vertical advection are complicated due to the various circulation patterns at these altitudes, for example, by the mesospheric residual circulation (e.g., Garcia & Solomon, 1985), or the lower thermospheric residual circulation at 95–110 km (e.g., Qian and Yue (2017), or both, depending on the season, altitude, and latitude. Further, Jones et al. (2018) demonstrated that the thermospheric spoon mechanism is most effective in the altitude regime where the thermosphere is transitioning from a fully mixed state to a diffusively separated state, indicating that the vertical advection/transport effects may also contribute to this uncertainty. Finally, as acknowledged above, the total vertical flux of atomic oxygen calculated for in Equations 1 and 2, is really the aggregate effect of the eddy, molecular, and “bulk” vertical wind on O, and thereby ultimately driving uncertainty in k .

However, modeling studies suggest that the coupling from above the MLT region is mainly due to molecular diffusion (see Figure 4b in Jones et al., 2017). For example, summing O densities above 140 km versus the columnar O density between 90 and 140 km leads to a ratio of ∼0.075 (or 7.5%). If we then assume, for example, a ∼30% change in O densities due to redistribution of O in the thermosphere, given diffusive equilibrium (dominated by molecular diffusion above 115 km), there would be less than a ∼2.3% change in the 90–140 km O density. That small change would have minimal effects on the determination of k .

In summary, the advection contamination at lower altitudes (below 90 km) where O lifetimes associated with recombination are short, vertical and horizontal advection contamination is minimal. At higher altitudes 90–96 km, advection can lead to uncertainties that approach 10%–20%, at the EQ. Given the discussion in the previous paragraph, uncertainties due to advection in the thermosphere would best be accomplished by refined global circulation models which would incorporate the details of the vertical distributions of temperature, composition, and dynamics. The contrast in the results between the EQ and NH/SH, where the strong SAO in the EQ zone and nearly absent in the NH, strongly suggest that even between 90‐96 km, contamination from horizontal advection effects between those two zones is minimal.

Additional Uncertainties and Future Studies

Figure 14 shows the near global (±55° average of k altitude profiles for SABER. The main difference between this study and S19) is an increase in the total profile by a factor of ∼2 at 83 km to 1.5 at 96 km. The primary reason for this is a lesser amount of O in the mean O profiles which leads to an increased k in order to meet the flux integral (Equation 2) from the 80 km lower boundary condition. Additional minor contributions are attributed to (a) the method of averaging k from the three zones, for the 26 periods in the annual cycle, rather than computing a global mean from the global mean O density profile, (b) MSIS2.0 that contains the relevant background atmosphere affecting the chemistry of O loss for each of the 78 temporal/spatial conditions, and (c) the O density was determined for nighttime conditions only, whereas the earlier study determined O density for day and night conditions. The O density uncertainties increase to a lower limit of 80 km, contributing to the uncertainty in k below 83 km (Figure 14). The chemical model was described by Mlynczak et al. (2013, 2018). Table 2 in Mlynczak et al. (2018) describes uncertainty of 20%, whereas the analysis by Mast et al. (2013) describes a smaller error. The uncertainty in k is less sensitive to O density but strongly sensitive to the gradient. A case study changing O density by 20% uniformly in altitude above 80 km resulted in a change in k of 1.6% at 96 km, increasing with decreasing altitude to 2.0% at 89 km, 3.0% at 84 km, and 6.3% at 80 km. The statistical uncertainties of 16 years of limb data inversions with geophysical variations of waves (gravity waves, tides, and planetary) contributions are unknown, but we believe less than those due to the uncertainties in the O density.

Clearly, the sophisticated approaches to define, track, filter, and propagate the meteorological effects from below as was done by Becker and von Savigny (2010); Grygalashvyly et al. (2012); A. Z. Liu et al. (2016), among others are evolving and constantly being improved. As the coupling of the atmosphere from the troposphere to the thermosphere is developed, the climatology of the minor constituents in the MLT, and parameterized transport effect from these climatologically driven studies, will play a role in refining and improving the process, with model validation.

Conclusions

The intra‐annual variation of eddy diffusion in the MLT region has been quantified using an analysis of the observed, average atomic oxygen profiles at 80–96 km by the TIMED‐SABER instrument (2002–2018) and 80–105 km by the Envisat SCIAMACHY instrument (2002–2014). The analysis method described for the determination of global mean k by S19 was used, along with the background atmospheric species were determined by the MSIS2.0 model (Emmert et al., 2021).

A list of our prominent findings are as follows:

The AO variation in k peaks in summer (near solstice) and is at a minimum in winter, with an amplitude factor change of ∼4X between solstices, in both the NH and SH at 83 km.

The difference between the winter and summer losses of O require an AO amplitude in the downward flux (nv) at 96 km, for 10%. The intra‐annual variation of the eddy diffusion velocity (v) is invariant at that altitude, but the O density is not, with measurements and analysis by Chen et al. (2019). The minimal O density below 90 km in the mid‐latitude summer is consistent with enhanced depletion and chemical O loss in summer.

The analysis of k in the EQ region resulted in a large SAO amplitude of 25% from solstice to spring equinox, and less from solstice to fall equinox.

The EQ k annual‐average profiles are uniform with altitude, with a value of 1.1 × 106 cm2s−1. The vertical extent of the large k , the upper limit of the data at 96 km, likely contributes to the simultaneous increase in O density, rather than the opposite effect, observed at mid‐latitude with the AO cycle. The uniformity with altitude and significantly larger SAO amplitude support the hypothesis that wave‐tide coupling contributes to turbulence and wave mixing at all EQ altitudes. The stratospheric wind minimum at the EQ also likely minimizes the filtering of waves reaching the MLT.

The climatology of the MLT k supports the Qian et al. (2009) annual mean, but not the intra‐annual variability in the downward flux at 97 km, similar to what was reported by Salinas et al. (2016). The variation in the AO (at solstice) is partially due to the variability in the O density, and consequently, the downward flux, which is highly variable with season. The intra‐annual global mean k and O density in global models would be best replaced with seasonal effects for the respective hemisphere, in order to better represent coupling effects taking place in the MLT.