Thermospheric neutral winds present a highly dynamic behavior with changing geophysical conditions. Understanding their variability becomes critical, as they transfer energy and momentum in the upper atmosphere and directly or indirectly affect the dynamics, morphology and composition of the ionosphere, which can disrupt radiocommunication and navigation systems (e.g., Wang et al., 2021).

Waves present in the thermosphere play a particular role and add significant variability in the wind system. Generally, these waves are classified as (i) planetary waves, which are global in scale with periods of up to several days to a month (Forbes, 1996); (ii) tides, which are also global in scale but with periods that are sub-harmonics of solar and lunar days (Oberheide et al., 2015); and (iii) gravity waves, which are medium-to small-scale oscillations with periods ranging from a few minutes to several hours (Fritts and Alexander, 2003).

Migrating and non-migrating tides originating from the lower atmosphere are recognized as important players in the vertical coupling between the lower and upper atmosphere. The general morphology of these tides has been studied extensively over the past decades (e.g., Lindzen, 1981; Teitelbaum and Vial, 1981; Miyahara et al., 1993; Liu et al., 2010; Forbes et al., 2017) and can be observed as global oscillations in winds, density, temperature and other atmospheric fields. In general, they transfer energy and momentum from the lower regions of the atmosphere into the upper atmosphere and generate longitudinal variations in the various atmospheric state parameters and can modify the global ionosphere-thermosphere (I-T) system (e.g., Immel et al., 2006; Forbes, 2007). Therefore, it is important to understand the spatial and temporal variability they generate in the upper atmosphere.

Thermospheric data derived from satellite accelerometers (e.g., Bruinsma and Biancale, 2003; Sutton et al., 2007; Doornbos et al., 2010) has been used to study global longitudinal structures produced by non-migrating tides. For example, Forbes et al. (2012) used densities from the SETA, CHAMP and GRACE satellites to investigate vertical tidal propagation and to identify the tidal oscillations in the longitudinal structures present in the data; Häusler and Lühr (2009) investigated the non-migrating tidal spectra in zonal wind data at equatorial latitudes obtained from the CHAMP accelerometer with an emphasis on the annual variation of the wave-4 structure at 400 km altitude; Lieberman et al. (2013a) investigated tidal variations in longitudinally averaged CHAMP global zonal winds; Gasperini et al. (2015) used GOCE neutral densities and zonal winds and TIMED-SABER temperatures to study vertical coupling in the thermosphere; Liu et al. (2016) found the presence of wind jets aligned with the magnetic equator in GOCE zonal winds; Dhadly et al. (2020) studied the latitudinal variation in intra-annual oscillations in the GOCE cross-track neutral winds.

This study focuses on longitudinal structures present in zonal winds produced by non-migrating atmospheric tides and observed by the GOCE satellite during dawn and dusk. The year-to-year progression during December solstice is investigated and the contributions from zonal wave-1 to wave-5 structures is studied. The GOCE results are compared to the Climatological Tidal Model of the Thermosphere (CTMT) (Oberheide et al., 2011).

Throughout this paper, the standard nomenclature for tides is used; DEs (DWs) is an eastward (westward) propagating diurnal tide with zonal wavenumber s. For semidiurnal tides, an S is used in place of the D.

This paper is organized in the following manner. Section 2 describes the data used in the study, while Section 3 describes the methodology. In Section 4 the combined contributions of zonal wave-1 to wave-5 structures are presented, and their individual contributions are shown in Section 5. In Section 6, the GOCE results are compared to the CTMT model. Finally, Section 7 provides a summary and discussion of the results.

The Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite was launched on 17 March 2009 into a dawn-dusk Sun-synchronous orbit with an inclination of 96.7° (near-polar orbit) and an altitude of ∼260 km. Its main objective was to study Earth’s gravity field, with thermospheric densities and winds calculated later combining accelerometer and ion thruster data, together with GPS tracking and star camera data. For a description of the determination algorithm see Doornbos et al., 2010. The data set version 2.0 was used for this study, which had been reprocessed with a new implementation of the algorithms (Visser et al., 2019). The algorithm uses a new satellite geometry and aerodynamic model representation (March et al., 2019a), with a new setting of the aerodynamic energy accommodation coefficient (March et al., 2019b).

Even though the GOCE satellite altitude is on the average ∼260 km, it varies with time, starting initially at ∼270 km in 2009 and decreasing to ∼250 km in 2013. The local time corresponding to dawn and dusk also varies. In 2009 the local solar time at the equator crossing for dusk is ∼18 h and by the end of the mission it reaches ∼19 h.

The errors in the GOCE zonal wind are of the order of ∼10%–20%, with the dominant source of errors being biases due to instrument calibration and external models used in the calculation of the winds (Doornbos et al., 2010). In this work, the errors will be attenuated by using residuals, which account for these biases.

Due to the availability of geomagnetically quiet-time data, the December solstice was selected to study the year-to-year progression of the longitudinal variability in the GOCE zonal winds. In order to calculate the perturbations in the GOCE zonal wind measurements, a 27-day window of data was selected for each year; this helps to minimize the effect of the Sun’s rotation. These selected windows are centered as close as possible to the December solstice for the years 2009, 2010, 2011, and 2012, taking into account gaps in the data and geomagnetically active periods. For reference, figures for June 2010 and 2011 are also included in the Supplementary Material. The selected windows are:

• December solstice 2009: 12 December 2009 to 07 January 2010

Figure 1A shows the daily F10.7 cm radio flux for the duration of the GOCE wind data set. The periods selected are indicated by green (December solstices) and purple (June solstices) vertical bars. The average F10.7 value for the selected window in December 2009 was 76 sfu. For 2010 that value was 81 sfu. The averages for 2011 and 2012 were higher, at 132 sfu and 109 sfu respectively. The mean F10.7 value for June 2010 was 75 sfu and for June 2011 it was 95 sfu. The geomagnetic activity in these periods was, in general, very low. The average Kp values did not exceed 1⁰ during the December solstice windows and it was below 2⁻ for the June periods. With the exception of one day in December 2011, the daily Kp index is always below a value of 3 for all of the December periods, and reaches 3.3 once during each of the selected June periods. The reason for not including June 2012 and June 2013 in our analysis is due to data gaps and the presence of geomagnetic storms during these periods which precluded us from finding suitable 27-day windows. Figure 1B shows, for the same period as above, the Oceanic Niño Index (ONI), a 3-month running mean of sea-surface temperature anomalies (Bamston et al., 1997) with respect to the mean from 1971 to 2000 in the region 120 to 170°W and 5°N to 5°S. The ONI is used to classify El Niño Southern Oscillation (ENSO). El Niño conditions are present when ONI exceeds +0.5K for 5 consecutive months whereas La Niña conditions correspond to values of ONI of −0.5K or lower. ENSO is categorized into weak (ONI values of 0.5–0.9), moderate (ONI 1.0-1.4), strong (ONI 1.5-1.9) and extreme (values of ONI 2 or higher). Figure 1C shows the Quasi Biennial Oscillation (QBO) U30 index, which corresponds to the zonally averaged wind at 30 hPa over the Equator. It is observed that throughout the span of the GOCE thermospheric data set both of these indices exhibit strong variations which will be discussed in Section 7.

Following a similar approach as outlined in Molina (2022) as well as in the companion paper by Molina and Scherliess (2023) the data corresponding to each 27-day window were separately analyzed for dawn (∼06 h local time) and dusk (∼18 h local time) conditions and sorted into bands of 1° of geographic latitude from −50° to 50° N. In each latitude band the median value was calculated. Next, the zonal wind perturbations (to be called dWind) were obtained for every GOCE zonal wind observation by computing: d W i n d = W i n d − M e d i a n , where the median subtracted from each wind measurement is the value that corresponds to the latitude band where the measurement is located.

For the longitudinal analysis the spatial resolution of the data was further reduced by organizing the dWind values into bins of 5° of latitude and 2° of longitude. In each bin, data points that fell beyond two standard deviations from the corresponding mean were filtered out and the remaining data points were averaged. Figure 2 shows, as an example, a global map of the calculated dWind for the December 2009 period. Figures 2A, B show each individual dWind data point whereas Figures 2C, D show the binned and averaged data. It can be seen that our binning and averaging preserves the main characteristics of the perturbations but attenuates the noise in the data. A longitudinal structure can be recognized from the plots as well, both in the dawn and the dusk perturbations. This structure presents alternating bands of positive (eastward) and negative (westward) wind perturbations that are the main focus of this paper.

In order to analyze the longitudinal structures and to elucidate their underlying wave characteristics, a 1-D Fast Fourier Transform (FFT) was applied to the binned and averaged deviations separately in each latitude band. Here, the 1-D FFT will provide the different longitudinal frequencies present in the zonal wind perturbations. Because the 1-D FFT is applied separately in each 5° latitude band, the results for each band will be independent from each other.

For a fixed local solar time, the tidal perturbations T ∼ can be expressed as (Oberheide et al., 2003): T ∼ = ∑ s , n T s , n ⁡ cos ω n t − t s , n − λ s + n where T s , n is the amplitude, ω n is the wave frequency, t is the local solar time, t s , n is the time of maximum amplitude with respect to 0° longitude, λ is the longitude, s is the zonal wavenumber and n is the number of cycles per day ( n = 1 for the diurnal components, n = 2 for the semidiurnal). This equation implies that when observing at a constant local time, the observed zonal wavenumber corresponds to s ′ = s + n . Depending on the sign of s , the wave propagation is eastward for positive values of s and westward for negative values of s .

Since the GOCE data are analyzed separately for dawn and dusk and thus correspond approximately to a fixed local time of ∼06 h for dawn and ∼18 h for dusk, the longitudinal variations in the GOCE data will appear as waves where the contributions from individual tides cannot be separated. Therefore, for example, both SE2 and DE3 will appear as part of a zonal wave-4 structure in the data, and their contributions will be combined in the FFT results. Since for migrating tides = − n , and consequently s ′ = 0 , their contribution appears as a constant in the tidal perturbations and becomes part of the median that we had subtracted from the zonal winds in our analysis. Consequently, the contribution of migrating tides is largely eliminated, and our results pertain to only non-migrating tides.

Figure 3 shows an example of the amplitude spectrum obtained from the FFT for the December 2011 period. Shown are the wave-m components up to wave-9. As already noted, these wave-m components do not pertain to a certain tidal wave, but instead are a combination of multiple tides. It is evident from Figure 3 that in this case the largest amplitudes are found in the first three wave components.

Figure 4 shows the wave amplitudes for December dawn conditions separately for each year from 2009 to 2012 (Figures 4A–D). Shown are the amplitudes for wave-1 to wave-9 as a global average (blue bars), averaged over the northern hemisphere (red bars), and averaged over the southern hemisphere (yellow bars). Here, the global values were obtained by averaging the individual wave-amplitudes in each latitude bin from −45° to 45° geographic latitude, and the northern/southern hemisphere averages were obtained by averaging the corresponding values from 0° to 45° and from 0° to −45°, respectively. It is interesting to note that during December 2009 the wave-1 amplitude in the northern hemisphere is more than twice the value found in the southern hemisphere. A similar hemispheric asymmetry can also be seen for wave-2 and wave-3 during December 2012. Figure 5 shows the same as Figure 4, but for dusk conditions. During this time, hemispheric asymmetries are present for wave-3 for December 2011 and for wave-2 and wave-3 during December 2012. Otherwise, the northern and southern averaged amplitudes are comparable to each other.

Figures 4, 5 also show that most of the wave amplitudes are concentrated in the first few wave numbers. In fact, our analysis shows that 75%–85% of the variability is represented by the first five components. As a consequence, we have limited our further analysis to only consider these first five wave-numbers.

Specifically, we have initially investigated the combined effect of wave-1 to wave-5 applying an inverse FFT (IFFT) after filtering out all contributions with zonal wavenumbers higher than 5. This was followed by a study of the individual effect of each wave-m component up to wave-5. For this all components of the FFT except the one corresponding to that particular wavenumber were filtered before the IFFT was performed. This step was repeated for all wave components from wave-1 through wave-5. In the following the combined contributions from wave-1 to wave-5 will be presented followed by a presentation of the individual contributions from wave-1 to wave-5.

The total December GOCE zonal wind perturbations (dWind) obtained from the IFFT produced by the contributions from wave-1 through wave-5 are separately shown for each year from 2009 to 2012 as a function of latitude and longitude in the first four rows of Figure 6. Figures 6A–D correspond to dawn and Figures 6E–H to dusk conditions. Positive (negative) values are shown as red (blue) colors and correspond to eastward (westward) perturbations. Figures 6I, J show the corresponding model results obtained from CTMT that will be discussed in Section 6.

In general, the total GOCE zonal wind perturbations at dawn remarkably resemble each other from one year to the next. During dusk, the GOCE zonal wind perturbations also show a good agreement for the low solar flux years 2009 and 2010 and as well as for the higher solar flux years 2011 and 2012. However, a clear change in the global pattern can be seen from 2010 to 2011. In the following sections, the global perturbation patterns are presented and agreements and differences from year to year are described in more detail.

In the following sections the individual contributions of wave-1 to wave-5 to the total zonal wind perturbations are analyzed separately. The results are first shown for dawn followed by the corresponding results for dusk.

In this section we will compare the GOCE results with simulations obtained from the Climatological Tidal Model of the Thermosphere (CTMT) (Oberheide et al., 2011) model. CTMT is an observation-based model that includes amplitudes and phases for six diurnal (DW2, DW1, D0, DE1, DE2, DE3) and eight semidiurnal (SW4, SW3, SW2, SW1, S0, SE1, SE2, SE3) tidal components for temperature, density, zonal, meridional and vertical winds from 80 to 400 km of altitude, pole-to-pole, and for moderate (F10.7 = 110 sfu) solar flux conditions. The model is based on Hough Mode Extensions to mean tidal diagnostics obtained from the TIMED Doppler Interferometer (TIDI) and the Sounding the Atmosphere using Broadband Emission Radiometry (SABER) instruments onboard the Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics (TIMED) satellite. CTMT perturbations have been compared to satellite observations (Forbes et al., 2012; Lieberman et al., 2013b; Forbes et al., 2014; Forbes et al., 2022), used as boundary conditions for numerical experiments (e.g., Jones Jr et al., 2019) and to interpret and explain ground-based observations (e.g., Yuan et al., 2014).

The monthly CTMT amplitudes and phases are provided on a latitude/altitude grid with 5° latitude and 2.5 km altitude resolution. For our comparison of the GOCE zonal wind perturbations with those predicted by CTMT we have calculated for each individual GOCE data point the corresponding CTMT value. Specifically, we have interpolated the provided amplitude and phase data to the same days and locations of each individual GOCE observation using a linear interpolation in each dimension. Once the corresponding CTMT amplitudes and phases were obtained, the corresponding CTMT zonal wind perturbation for each individual tide was calculated using: T ∼ n , s = A n , s ⁡ cos n Ω t L − s + n λ − ϕ n , s where T ∼ n , s is the tidal perturbation, n is the subharmonic of a solar day (n = 1 for diurnal tides and n = 2 for semidiurnal tides), s is the zonal wavenumber, A n , s is the amplitude, Ω is the rotation rate of the Earth, t L is the local solar time, λ is the longitude and ϕ n , s is the phase. Finally, the individual tidal components were added to obtain the total simulated zonal wind perturbation.

When investigating the longitudinal variations in the CTMT simulated zonal wind perturbations at a fixed local time the contributions from migrating tides appear as a constant and longitudinal variations are only due to non-migrating tidal components that will appear as a wave-m structures spanning from wave-1 to wave-5. In particular, D0, DW2, SW1 and SW3 will appear as a wave-1. Wave-2 will be comprised of DE1, S0 and SW4. Wave-3 will be composed of DE2 and SE1. DE3 and SE2 will appear as a wave-4. SE3 is the only component in CTMT that will appear as a wave-5.

Once the CTMT perturbations were calculated, the same steps performed on the GOCE zonal wind perturbations described in Section 3 were applied to the simulated CTMT zonal wind perturbations.

We have used low-to mid-latitude zonal wind observations from 2009 to 2012 obtained by the GOCE satellite near 260 km altitude during geomagnetically quiet times to investigate the interannual variation of the longitudinal variability of the zonal wind near dawn and dusk. The focus of the study was to investigate the year-to-year progression of the longitudinal variability in the GOCE zonal winds during December solstice produced by nonmigrating atmospheric tides. For each year the contributions from zonal wave-1 to wave-5 structures were separately determined and compared. The GOCE longitudinal zonal wind variations were also compared with corresponding results obtained from the CTMT model.

To determine the longitudinal variation in the GOCE zonal wind measurements, a 27-day window (one solar rotation) was selected for each year and the data were separately analyzed for dawn and dusk from −50° to 50° geographic latitude. The temporal and spatial variability of the zonal winds was then presented as a variation about the zonal mean values and decomposed into its underlying wave-m structures using a Fourier analysis. This approach largely eliminated the contribution of migrating tides from our analysis but also precluded us from separating the longitudinal variations into their individual tidal components. Consequently, our zonal wind perturbations appear as zonal wave-m structures that result from the superposition of individual non-migrating tidal components. It is important to highlight that in this methodology, the Fourier analysis is applied to each latitude band individually; therefore, the results for each one of these bands are independent of each other. The coherence between the adjacent latitude bands presented in Figures 6–8 gives confidence in the structures observed.

It was found that 75%–85% of the longitudinal zonal wind variability could be explained as due to waves that generate zonal wave-m structures with m up to 5, and therefore, this was the cutoff used in our Fourier analysis. A clear interannual progression of the individual wave components could be observed in the resulting structures.

In general, the total GOCE zonal wind perturbations at dawn remarkably resemble each other from one year to the next. During dusk, the GOCE zonal wind perturbations also show a good agreement for the low solar flux years 2009 and 2010 as well as for the higher solar flux years 2011 and 2012. However, a clear change in the global pattern can be seen from 2010 to 2011. Note that the F10.7 solar flux values during 2009 and 2010 were lower than during 2011 and 2012 and the change in the pattern might be related to the change in solar flux. A similar good agreement in the total GOCE zonal wind perturbations can also be observed for June 2010 and June 2011 (shown in Supplementary Figure S1). As pointed out above, the average F10.7 solar flux during these two June periods also only differs by 20 sfu.

Both ENSO and the QBO are also known to produce year-to-year variability on atmospheric tides (e.g., Oberheide et al., 2009; Warner and Oberheide, 2014). As mentioned above, the GOCE results for each December pair 2009-2010 and 2011-2012 are similar. Figure 1B shows, however, that the ONI values for 2009 and 2010 are very different and correspond to moderate El Niño and strong La Niña conditions, respectively. December 2011 presents moderate La Niña conditions, and December 2012 does not present either. From Figure 1C, it is also evident that the QBO U30 index is different for the December 2009 and 2010 pair. Differences are also observed in the index for the June 2010 and 2011 pair. Based on these observations, it appears that the year-to-year variations seen in GOCE are not the result of either variations in ENSO or the QBO. During dawn, the total December solstice GOCE zonal wind perturbations at low- and mid-latitudes range from about −40 m/s to 50 m/s and exhibit a clear longitudinal structure consisting of four distinct bands. These bands alternate between eastward and westward perturbations and are generally tilted westward with increasing latitude (north-west alignment). The individual wave-m structures during dawn are also remarkably similar from year to year in both phase and magnitude, but also exhibit differences when comparing the individual wind perturbations for the low solar flux years (2009-2010) to those of the higher solar flux years (2011-2012).

During dusk, the low and mid-latitude GOCE zonal wind perturbations for December solstice also range from −40 m/s to 45 m/s. Like dawn, a clear longitudinal banded structure can be identified at with four bands alternating between eastward and westward perturbations. This time, however, the bands are generally tilted eastward with increasing latitude (north-east alignment). Similar to the general morphology during dawn, the individual wave-m structures are again similar from year to year in both phase and magnitude, but this time the differences between the individual wind perturbations for the low solar flux years (2009-2010) and those for the higher solar flux years (2011-2012) become even more evident. Here again, a similar result is found for the two June periods where the individual wave-m structures resemble each other year-to-year (see Supplementary Figure S2 for dawn and Supplementary Figure S3 for dusk).

Some of the characteristics just described are observed in a companion paper by Molina and Scherliess (2023), who have used the same GOCE zonal wind observations to determine their spatial and temporal correlations. In particular, it is mentioned that the longitudinal/temporal correlations for December 2009 indicate that the structures that generate them have a north-west orientation during dawn (Figure 10 in their paper) and north-east orientation during dusk (their Figure 11). It is noteworthy that this is observed in Figure 6 of this paper as the general orientation of the longitudinal bands for December 2009 for dawn and dusk, respectively. The zonal wave structures up to wave-5 obtained in this work were also subtracted from the associated perturbations, and the longitude/time correlation coefficients were subsequently calculated for their study. The resulting correlation coefficients (Figures 12 to 15 in their paper) suggest that the zonal wave-m structures described in this paper are largely responsible for the original patterns observed in their longitudinal/temporal correlations.

A comparison of the GOCE results with simulations obtained from the CTMT model has revealed a very good agreement between the CTMT and the December 2009 and 2010 GOCE zonal wind perturbations at dawn, especially in the northern hemisphere, where the observed longitudinal bands align within 10° of those seen in the model results. However, even though there are striking similarities between the general morphology of the GOCE and CTMT zonal wind perturbations, there are also noticeable differences. During dusk the CTMT zonal wind perturbations show significant differences with the observed GOCE zonal wind perturbations. The structures present in the CTMT perturbations do not resemble the general morphology of the GOCE results for any of the years. In particular, the CTMT wave-1 component significantly differs from the GOCE results (this can also be seen in Supplementary Figures S2, S3 for June solstice). There could be several reasons for this discrepancy. Oberheide et al. (2011) compared the CTMT densities near 400 km altitude with those obtained from the CHAMP satellite and found that the agreement with the tidal components DW2, D0, SW1, and SW3 was poor. As mentioned before, these are also the components that constitute the wave-1 structure. They suggest that the observed differences are the result of the generation mechanisms of these tidal components, which are believed to be generated by hydromagnetic coupling between various waves (Jones et al., 2013) which is not captured by the CTMT formulation. Part of the discrepancy between GOCE and CTMT at dusk might also be the result of our use of a geographic coordinate system in our analysis. Liu et al. (2016) found the presence of wind jets aligned with the magnetic equator in GOCE zonal winds. Zhang et al. (2018) used CHAMP zonal winds together with TIEGCM model simulations to investigate the effect of the geomagnetic field on longitudinal variations observed in zonal winds. They conclude that the large-scale longitudinal variations are produced by the geomagnetic field structure and might be the result of temporal variations of ion drag and pressure gradient forces. Indeed, the variation captured by TIEGCM in their Figure 10 is reminiscent of a wave-1 structure for both dawn and dusk.

Regarding the differences in magnitudes of the CTMT wind perturbations, which are only about one-third of the corresponding GOCE values, it needs to be noted that it is not clear whether this difference is the result of an underestimation of zonal wind perturbations by CTMT or due to a systematic overestimation of the GOCE wind data or both. Jiang et al. (2021) for example, compared GOCE zonal winds with wind observations obtained from ground-based Fabry-Perot Interferometer (FPI) measurements at low and mid latitudes and reported an overall overestimation of the GOCE winds when compared to the ground-based data. They found that the magnitudes from GOCE are generally larger than the FPI winds by a factor of 1.37–1.69, consistent with FPI comparisons at high latitudes where factors of 1.2–2.0 were reported by Kärräng (2015), albeit using an earlier version of the GOCE data. An earlier version of the data was also used by Dhadly et al. (2017, Dhadly et al., 2019) who reported a magnetic latitude-depended bias in the GOCE data when compared to WINDII, SDI and FPI observations at high latitudes. The possible presence of a bias in the GOCE data, however, would have only affected our zonal mean values and consequently would have been eliminated when calculating the wind perturbations. A factor difference, however, would also affect our perturbation values and consequently our reported values should be scaled by this factor.

Finally, it should be noted that our analysis only pertains to December and June solstice conditions at dawn and dusk and an investigation of the year-to-year progression of the longitudinal variability of the zonal wind during other seasons and local times needs to be performed in the future.