Quiet Time Thermospheric Gravity Waves Observed by GOCE and CHAMP

The Gravity Field and Steady‐State Ocean Circulation Explorer (GOCE) and CHAllenging Minisatellite Payload (CHAMP) satellites measure in‐situ thermospheric density and cross‐track wind. When propagating obliquely to the satellite track in a horizontal plane (i.e., not purely along‐track or cross‐track), gravity waves (GWs) can be observed both in the density and cross‐track wind perturbations. We employ the Wavelet Analysis, red noise model, dissipative dispersion and polarization relations for thermospheric GWs, and specific criteria to determine whether a quiet‐time (Kp < 3) thermospheric traveling atmospheric disturbances (TADs) event is a GW or not. The first global morphology of thermospheric GWs instead of TADs is reported. The fast intrinsic horizontal phase speed (cIH > 600 m/s) of most GWs suggests that they are not generated in the lower/middle atmosphere (where cIH < 300 m/s). A second population of GWs with slower speeds (cIH = 50–250 m/s) in GOCE are likely from the lower/middle atmosphere, but they occur much less frequently in CHAMP. GW hotspots occur during the high‐latitude and the winter midlatitude regions. GW amplitudes exhibit semi‐annual and annual variations. These findings suggest that most GOCE and CHAMP GWs are higher‐order GWs from primary GW sources in the lower/middle atmosphere. Finally, the average propagation direction of the CHAMP GWs exhibits a clear diurnal cycle, with clockwise (counterclockwise) occurring in the northern (southern) hemisphere and equatorward propagation occurring at ∼13 LST. This suggests that the predominant GW propagation direction is opposite to the background wind direction.

XU ET AL.
The influence of GWs becomes pronounced in the middle and upper atmosphere, impacting the variability, winds and temperatures, and global circulation.It is well known that GWs serve as drivers for global-scale phenomena such as the quasi-biennial oscillation (Baldwin et al., 2001) and the mesospheric residual circulation (Holton, 1982).Recently, there has been a growing interest in understanding the connection between lower and middle atmospheric activities, such as the effects that GWs from the polar vortices have on the variability of the thermosphere and ionosphere.For instance, studies have revealed that during periods of weakened polar vortex states, the upper stratosphere and lower mesosphere experience reduced GW activities (Harvey et al., 2023), resulting in weaker thermospheric GW activities (Becker, Vadas, et al., 2022) and fewer traveling ionospheric disturbances (TIDs) (Frissell et al., 2016).These findings underscore the intricate interplay between lower atmospheric dynamics and the behavior of GWs throughout the atmosphere, offering insights into the broader impacts of atmospheric processes from below.
GWs are also important to telecommunication and satellite operation.Perturbations in the thermosphere (e.g., the occurrence of equatorial spread F, also called equatorial plasma bubbles) have received intensive attention because of their complex physics and severe influence on telecommunication and navigation (Aa et al., 2020).The seeding of those phenomena is the least understood, and GWs have been commonly suggested to play a key role in the seeding process (H.X. Liu et al., 2017;Takahashi et al., 2009).
Although thermospheric GWs play a vital role in thermospheric dynamics, their measurement in this region is challenging.Previous studies have employed indirect methods to study thermospheric GWs.For example, Bhatt et al. (2023) utilized airglow perturbations captured by optical sensors like all-sky imagers and Fabry-Perot Interferometers to derive nighttime wind and temperature perturbations.Another approach uses Global Positioning System (GPS) Total Electron Content data sets to retrieve GW characteristics from TIDs around the F2 layer (Azeem et al., 2015;Nishioka et al., 2013;Xu et al., 2019).These data have high horizontal and temporal resolution, but are limited to regions with dense ground-based receiver networks, such as the continental US.Satellite-based remote sensing methods have also provided valuable observations in the thermosphere/ionosphere; however, they are subject to limitations of low time or horizontal resolution.For instance, the TIMED Doppler Interferometer (TIDI) onboard the Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics (TIMED) satellite and Michelson Interferometer for Global High-resolution Thermospheric Imaging (MIGHTI) onboard Ionospheric Connection Explorer (ICON) satellite measure horizontal wind profiles along the satellite orbit but with coarse horizontal resolution (∼750 km for TIDI (Yee et al., 2003) and 250-500 km for MIGHTI (Cullens et al., 2022;Dhadly et al., 2021)); The Global-scale Observations of the Limb and Disk (GOLD) mission can observe neutral temperature disturbances at an altitude of ∼150 km, but with coarse horizontal (250 × 250 km 2 at nadir) and time resolutions (∼30-min cadence for day disk images) (McClintock et al., 2020).To overcome these limitations, in-situ spacecraft missions equipped with precise accelerometers have emerged as powerful tools for studying traveling atmospheric disturbances (TADs) and thermospheric GWs.Low Earth orbit (LEO) missions such as CHAMP, GOCE, and Gravity Recovery and Climate Experiment (GRACE) have provided unprecedented accuracy and resolution in density and cross-track wind measurements (Doornbos, 2012).Previous studies have detected TADs, which include GWs, using accelerometer data onboard several spacecrafts (e.g., Bruinsma & Forbes, 2008;Li et al., 2023;H. Liu et al., 2017;Park et al., 2014Park et al., , 2023;;Trinh et al., 2018).However, due to the absence of along-track wind measurements, these studies have primarily focused on either density or cross-track wind disturbances along the satellite track.Whether these TADs are GWs or not is not known without detailed investigation of each event.
The paper is arranged as follows: In Section 2, we describe the GOCE and CHAMP data sets.In Section 3, we introduce (a) the methodology to pick out TADs through rectified Wavelet analysis and (b) the methodology to determine the GW intrinsic parameters (λ H , λ z , c H , τ Ir and their direction of propagation) using the dissipative GW dispersion and polarization relations, and to verify those TADs that are GWs by these solved intrinsic parameters.In Section 4, we show statistical (seasonal/monthly) results of thermospheric GWs obtained from GOCE and CHAMP data sets.Section 5 provides conclusions and an outlook.

GOCE and CHAMP Data Sets
The GOCE satellite was a European Space Agency mission that aimed to improve our understanding of Earth's gravity field and geodynamics.Launched on 17 March 2009, GOCE orbited Earth at an unusually low altitude of ∼250-290 km with an inclination of ∼96.7°, and local solar time (LST) of its ascending node drifted slowly from 18:13 LST (2009 November) to 19:33 LST (2013 October) (Doornbos, 2019).To maintain its low altitude for four to 5 years, the satellite used an ion thruster assembly that aligned with the orbit direction.GOCE carried a payload Electrostatic Gravity Gradiometer that consisted of six highly accurate accelerometers with a noise level of 10 −12 m s −2 Hz −1/2 (Doornbos et al., 2009;Visser et al., 2018).The orbital period of GOCE was ∼90 min.
The CHAMP mission was led by the GeoForschungsZentrum (GFZ, or Geo-research Centre) and launched on 15 July 2000, with the primary objective of measuring Earth's gravity and magnetic field to high precision.CHAMP was in a near-circular orbit with an inclination of ∼87°, and its ascending and descending nodes covered all LST (24-hr) about every four months.The satellite's altitude decreased gradually from 450 km at the beginning of the mission to 260 km in August 2010.The orbital period of CHAMP was about ∼90 min, with a slight decrease of ∼5 min due to its descending altitude during its life span.Equipped with a Spatial Triaxial Accelerometer for Research (STAR) accelerometer, the measurements of CHAMP were less precise than those of GOCE.The STAR accelerometer had a sensor in the axis parallel to the Z-axis of the satellite body-fixed XU ET AL. 10.1029/2023JA032078 4 of 35 frame (nominally pointed toward the satellite's nadir) that suffered from a malfunction, and was less sensitive than the sensors along the other two axes, but STAR still had a resolution <3 × 10 −9 m s −2 Hz −1/2 (Rodrigues et al., 2003), which was adequate for deriving in-situ density and cross-track wind measurements and TAD analysis.For both GOCE and CHAMP, we specifically select TAD and GW events that occurred when both concurrent Kp and the average Kp over the past 6 hr were less than 3, aiming to minimize wave events excited by geomagnetic disturbances.Doornbos et al. (2010) described the methodology for deriving in-situ density and cross-track wind measurements for the CHAMP mission.The approach for retrieving the density and cross-track wind for GOCE were essentially identical (Doornbos et al., 2013).First, when initially assuming a co-rotating thermosphere with density but no wind suggested by the Mass Spectrometer and Incoherent Scatter radar (MSIS) model, the modeled acceleration experienced by a spacecraft does not align with the measured acceleration.The discrepancy between the measured and modeled accelerations is mainly due to two parts: the directional discrepancy is attributed to the absence of cross-track wind in the model, the magnitude discrepancy is attributed to incorrect model density.To address this, a cross-track wind is introduced to the relative velocity (the difference between the bulk velocity of gas particles and the spacecraft velocity), and the absolute in-situ mass density is adjusted at each iterative step.This gradual iteration allows the modeled acceleration to converge toward the observed aerodynamic accelerations.It should be noted that the in-track wind is assumed according to a model value (Doornbos, 2019).However, given that the magnitude of in-track wind (∼ a few hundreds of m/s at most) is significantly smaller than the spacecraft velocity (∼8 km/s), and the low frontal area-to-mass ratio of the spacecraft, it was impossible to distinguish how much in total in-track acceleration discrepancy was induced by the in-track wind deviation and density deviation (Doornbos, 2012;Doornbos et al., 2010).In principle, the wind component in the vertical direction could be derived in a similar way as the cross-track wind, but there are two obstacles to overcome for CHAMP: (a) aerodynamic acceleration in the vertical direction is in general too small when compared with errors in the instrument calibration, the radiation pressure model, and the lift force model; (b) the malfunction of the CHAMP accelerometer in Z-axis prevents the acquisition of data in the vertical direction (Doornbos et al., 2010).Therefore, only the density perturbation and cross-track wind are available in the CHAMP data set analyzed in this paper.Because the GOCE mission carried more advanced accelerometers than CHAMP, its measurements have the capacity to derive both cross-track and vertical winds, for example, R. F. Garcia et al. (2014) introduced disturbances in both cross-track and vertical winds from the Tohoku earthquake and tsunami occurred in 2011.However, in the public-released the GOCE data set, the vertical component of the crosswind is assumed to be zero (Doornbos, 2019), the vertical wind perturbations are thus not publicly available in the GOCE and CHAMP data sets; instead, the density and cross-track wind data sets are available and are used in our study.
In this study, our primary focus lies on the density and cross-track wind observations obtained from GOCE during November 2009 to October 2013 and from CHAMP during 2004-2007, when satellites were at relatively stable altitudes (∼270 km for GOCE and ∼370 km for CHAMP).The data were sampled at a rate of 1 data point every 10 s (or 80 km in distance) for both GOCE and CHAMP, with the Kp index provided every 3 hr.Based on Doornbos et al. (2010), the rms density and wind errors in the data can be as low as 0.03% and 1 m/s, respectively.Although the systematic density error in CHAMP was estimated to be around ∼10%-15%, the same algorithm applied to GOCE data could potentially yield even smaller rms errors due to the presence of more and better instruments onboard.These favorable error characteristics allow us to study relatively small-amplitude and small-scale perturbations such as TADs and GWs.To provide accurate estimates of the errors associated with each observed wave event, we developed specific techniques, described in Section 3 of this paper.

Techniques for Diagnosing GWs From In-Situ Observations
The in-situ GOCE and CHAMP measurements provide valuable information about the amplitudes and scales of TADs along the satellite tracks.However, due to the presence of intrinsic noise in these measurements, it can be difficult to distinguish between noise and meaningful physical variations in the thermosphere.Despite this challenge, previous studies have successfully employed statistical methods to accurately represent the seasonal and global distribution of TADs from GOCE and CHAMP measurements (e.g., Forbes et al., 2016;H. X. Liu et al., 2017;Park et al., 2014;Trinh et al., 2018;Xu et al., 2021).In this study, we present a novel method to identify TADs from noise (Section 3.1) and then solve GWs from the TADs (Section 3.2).

Diagnosing TADs Using Rectified Wavelet Analysis Method and Red-Noise Model
Wavelet analysis is chosen as the preferred method in our study due to its ability to effectively handle localized and non-stationary wave packets, which are common in real-world scenarios.While the Fast Fourier transform (FFT) is suitable for quantifying perturbations with fixed wavelengths across short satellite tracks (e.g., Vadas et al., 2019), it may not work when wave packets have varying scales, frequencies, and phases, and when multiple wave packets are superimposed, making their separation challenging.In contrast, wavelet analysis offers a powerful approach to separate each wave packet based on its dominant period, amplitude, location, and duration.Wavelet analysis has been widely utilized by researchers, including Torrence and Compo (1998, hereafter TC98) and Chen (2016, hereafter C16).The versatility and effectiveness of wavelet analysis make it suitable for investigating the complex wave packets here.
The wavelet analysis derives a spectrum of variations arising from both noise and variations of TADs.To differentiate between perturbations with low confidence (noise) and high confidence (TADs), we assume that the noise can be modeled by red noise.The red noise model assumes that noise has a power (variance) spectral density that decreases as the scale decreases whereas the white noise model assumes that noise with different scales have equal power.In line with TC98, we first calculate the autocorrelation of the input series, which can be used to derive the expected normalized power of noise at each frequency.In this way, we can estimate the errors even though quantitative uncertainties for each data point are unavailable.Many geophysical time series can be modeled as either white noise or red noise (TC98).For instance, TC98 modeled the noises in the Niño3 sea surface temperature time series using the red-noise model.They then established significance levels and confidence intervals to diagnose wave packets with higher than 95% confidence level.
The univariate lag-1 autoregressive red noise model proposed by TC98 is a valuable technique for processing a large number of TADs, especially in scenarios where the input series is primarily governed by stochastic noise with different scales.Due to its effectiveness and computational efficiency under such assumptions, it serves as a practical and cost-effective method for automated processing.However, it is essential to recognize that the red noise model may not be suitable for all situations.In cases where the input series exhibits variations dominated by different scales of signals rather than noise, the red noise technique can lead to underestimating the confidence for TADs.Consequently, this may result in an underestimation of the number of TADs identified in the data set.Improving the estimation of confidence in each wave/noise packet remains a challenge.To address this limitation and achieve more accurate results, further investigations into the noise characteristics of the GOCE and CHAMP data sets are required.This additional effort will be critical for enhancing the robustness of our findings and ensuring the reliability of the outcomes obtained through the red noise model.
By following the practical guide provided in TC98, we calculate the wavelet power spectra  | ( track)| 2 for density and cross-track wind as functions of the data point index (n) and the along-track wavelength (λ track , which takes into account the satellite motion, see Appendix D for more details).These power spectra quantify the TADs along the satellite track.However, the TC98 method is biased toward larger scales when the input is considered as a superposition of monochromatic waves (Y.Liu et al., 2007;Maraun & Kurths, 2004).To reduce this bias, C16 replaced the wavelet function    () in TC98 with a modified wavelet function   mod =    .We adopted this change to obtain an unbiased wavelet spectrum W mod = 2C F W which we use to quantify the amplitudes and phases of TAD events in density and cross-track wind.Here, the correction factor C F is a function of λ track , and is defined in Appendix B.
In processing the GOCE and CHAMP data, it is common to encounter gaps in the measurements due to missing data.To address this issue and other computational considerations, we have divided the GOCE and CHAMP data sets into segments that are typically one full-circle orbit long.In Figure 1, we present the procedure for determining TADs in the density and cross-track wind data using a full-circle GOCE orbit on 5 July 2010 at 22:10-23:40 as an example.The full-circle GOCE orbit (red dots) begins in the westward direction (open red circle) and ends at the solid red circle (in the Northern Hemisphere).This orbit starts with the GOCE satellite in a descending mode and ends after completing the following ascending mode, which we refer to as the "DA (Descending-Ascending) mode."Due to the flattening of the reference ellipsoid representing the oblate Earth, the eccentricity of the orbit, and the perigee rotation rate, the altitudes of GOCE varies periodically (Doornbos et al., 2013).Figure 1b shows the in-situ density observation versus altitude within the track.From the relationship between the log-scaled density and height, we derive the thermospheric density scale height.To eliminate the density variation caused by the altitude change, we extrapolate the measured density XU ET AL.

10.1029/2023JA032078
6 of 35 to a fixed altitude, which is 270 km for GOCE and 370 km for CHAMP, based on the thermospheric density scale height in the observation.As shown in Figure 1c, the extrapolated density is denoted as ρ, which is the sum of the background density   (defined by the average of the density at the extrapolated altitude) plus the The green curve denotes the background density for wavelengths ≥3,000 km.(e) Unbiased amplitude spectrum of the detrended (λ track < 3,000 km) relative density perturbation (   ′ ∕ ).The thick black contour encloses TADs with >95% confidence for a red-noise process.Within the >95% confidence contours, peaks labeled with green crosses and green indices denote locations of TADs in density.Gray cross-hatched regions indicate the "cone of influence," where edge effects are important and waves are not trustworthy.(f-g) The same as (d-e) but for cross track wind perturbation   ′ track .In panel (g), peaks labeled with orange crosses and orange indices denote locations of TADs in cross-track wind.(h) Comparison between the sum of reconstructed perturbations    (blue, λ track < 3,000 km, confidence>95%) and the original detrended   ′ ∕ (black, λ track < 3,000 km).The horizontal axis only spans from 2 × 10 4 to 3 × 10 4 km.(i) The same as (g) but for reconstructed perturbations    (orange, λ track < 3,000 km, confidence>95%) and detrended   ′ track (black, λ track < 3,000 km).(j) The reconstructed perturbation pair of  10 * ρ (dark and light blue) and    (dark and light red) at λ track = 1,087.2km.The dark/light color denotes perturbations with confidence higher/lower than 95%.The "PS" means "phase shift" in units of degrees, "AR" means "amplitude ratio" in units of %*s/m.XU ET AL. 10.1029/2023JA032078 7 of 35 perturbation density ρ′. Figure 1d shows the relative density perturbation (black), defined as   ′ ∕ .Figure 1e shows the unbiased wavelet spectrum |W mod | in   ′ ∕ for TADs with λ track < 3,000 km in the relative density perturbation.The cross-hatched region shows the so-called "cone of influence" where edge effects are important.We only analyze TADs inside the center half of the orbit and outside of the "cone of influence".In the DA mode, TADs in the Southern Hemisphere are free from the edge effects on either end of the track, allowing us to focus on the peaks within the middle part of this circular orbit.To study TADs in the Northern Hemisphere, we re-process the same data set but divide the orbits by the "AD (Ascending-Descending) mode", so the "cone of influence" (endpoints) will be in the Southern Hemisphere.If a peak is located inside the thick black contours in Figure 1e, we have confidence higher than 95% that this peak is caused by a TAD rather than noise (Appendix B contains further information about the "confidence level of 95%".).For TADs in the cross-track wind, we directly extract the TADs from the observations since the wind is much less sensitive to changes in the altitude than the density.Figure 1f shows the cross-track wind perturbation (black), defined as   ′ track .Figure 1g B9) and (B12)).For each TAD event in our program, we apply a Monte Carlo simulation to generate an N-element complex array (N = 1,000 in current study) so that the N-element phasor array of    satisfies the 2D normal distribution in complex plane (a phasor is a complex number that represents a sinusoidal function.)The same process is applied to the phasor of   track to generate an independent stochastic N-element complex array.By dividing the complex array for    by the complex array for   track , we derive N PS values and N AR values (so-called "Monte-Carlo results" in this paper) between    and   track .In Section 3.2, we explain how we use these N pairs of PS and AR values to solve N possible GW solutions for each TAD event, taking the specific TAD event shown in Figure 1j as an example.Whether a TAD observed in the along-track data set is a GW or not (and if yes, how the GW parameters are distributed) will depend on the distribution of those N GW solutions.

Solving GWs Using Dissipative Dispersion and Polarization Relations of Thermospheric GWs
When a wave modulates the thermospheric density and wind simultaneously, phase shifts and amplitude ratios in the perturbations of two or more different parameters (e.g., u′, w′, ρ′, and T′) can be observed.Those observed perturbations can be used to determine parameters of the wave (e.g., λ H , λ z , ω Ir ).For example, Gross et al. (1984) estimated the directions of GW propagation within 180° by utilizing the phase shifts and amplitude ratios of O and N 2 perturbations in the data collected by the Atmospheric Explorer-C satellite.Innis and Conde (2002) utilized measurements from the Dynamics Explorer 2 satellite to determine GW propagation directions with an accuracy of within 90°; although their method included compressibility of air, it did not include viscous dissipation.In a later study, Vadas and Nicolls (2012) developed a generalized approach that accounted for full compressibility and realistic molecular viscosity in the GW dispersion and polarization relations, building on the original work by Vadas and Fritts (2005).In this section, we will introduce how to solve the intrinsic parameters of a GW event when it is measured as a TAD event in the density and wind by a satellite that flies through the thermosphere.
According to Vadas and Nicolls (2012), when the background conditions are known, both the sinusoidal perturbations in the relative density ( GW can be expressed as a function of horizontal and vertical wavelengths (λ H and λ z ) using the GW dispersion relation and polarization relations.We define their ratio (   ∕   ) as the complex function F: .Suppose we want to confirm whether TAD #2 identified in the density perturbation (see TAD #2 Figure 1e) is a GW event or not.The background parameters needed by the function F, such as N B (buoyancy frequency), c s (sound speed),   (density scale height), ν (kinematic viscosity,   = ∕ , where μ is the molecular viscosity), γ (ratio of specific heats), etc., are determined directly or indirectly from the Naval Research Laboratory MSIS (NRLMSIS) 2.1 model (Emmert et al., 2022).This model inputs the time and location for the TAD event.For TAD #2 event shown in Figure 1, (lat, lon) = (79.43°S,75.53°E) and UTC time = 2010-07-05_22:53:00.Figures 2a and 2b visualize the PS and AR of  ( ρ ρ ) for GWs with different (λ H , −λ z ) for the given background.The potential influence of uncertainties in the NRLMSIS2.1 model on GW solutions, as demonstrated by Vadas and Nicolls (2012), is acknowledged.However, the thorough quantification of this uncertainty necessitates a dedicated investigation, which lies outside the scope of the current study.
If GOCE and CHAMP had measured the horizontal wind perturbation    directly, then intrinsic parameters like λ H and λ z would have been directly determined by measurements of    and    .However, since GOCE and CHAMP only measured the cross-track wind perturbation   track (the perturbation perpendicular to the satellite track), λ H and λ z of the observed GWs cannot be determined directly unless the geometric relationship between the satellite and GW (i.e., the angle between the satellite path and the GW propagation direction in the horizontal plane) is known or determined.Figure 3 illustrates the geometry of a satellite moving at azimuth ψ, and a GW with  In Figure 3, λ ψ is the distance between the phase lines of GW along the azimuth ψ at an instant in time (λ H = λ ψ *|-cosϕ|), while the along-track wavelength λ track is determined from the peaks in wavelet spectra discussed in Section 3.1 (e.g., Figures 1e and 1g).If we assume that the satellite velocity is infinitely large, then λ ψ and λ track are identical.However, λ ψ and λ track differ somewhat due to the Doppler shift (especially for the fast GWs), and their relationship is derived in Appendix D, where U H is the background horizontal wind along the direction of GW propagation.Since λ track is obtained from the wavelet spectra and U H is derived from the HWM14 model (Drob et al., 2015), Equation 3 describes the relationship between λ H and ϕ. (Similar to NRLMSIS2.1, the current study does not contain a quantification of the uncertainty in the HWM14 model.Addressing this limitation is presently beyond the scope of our capabilities.)By plugging Equation 3and λ H = λ ψ * |cosϕ| into the function G, the number of unknowns on the right-hand side can be reduced from 3 (λ H , λ z and ϕ) to 2 (λ z and ϕ); thus the function G can be written as The complete form of function   is elaborated in Equation C11 of Appendix C.
The determination of the left-hand side of Equation 4 through the density and cross-track wind observations is explained in Section 3.1; however, solving ϕ and λ z is still a great challenge due to the computational complexity of linearizing function G within a reasonable processing time for automation purposes.An easier way to address this problem is to make two arrays for the right-hand side of Equation 4 and then utilize the arrays as a pair of "lookup tables" to find the best solution for ϕ and λ z with given PS and AR of  ( ρ ρtrack) .4a and 4b show the distribution of these solutions with perturbed PS and AR values; 86% of these solutions are in Sat.Q-II (ϕ∈(0 °, 90 °)).Therefore, this TAD event (observed by GOCE) is likely a GW event with a position in ϕ−λ z space in Sat.Q-II (ϕ∈(0 °, 90 °)).Since the satellite movement azimuth is ψ = −140.7 ° and the peak of ϕ = 43.6°(Figure 4c), the azimuth of this GW's propagation direction is ξ = ψ−ϕ = 175.7°. Figure 4d shows the location (green cross) and simultaneous propagation direction (magenta arrow) of the observed GW overplotted in a map.
Figures 1-4 illustrate the procedure for extracting TADs from along-track GOCE observations and determining if a TAD event represents a GW by solving its most probable intrinsic parameters.In Figures 4a and 4b, the concentration of Monte-Carlo results within Sat.Q-I provides strong evidence that this TAD is a GW.However, a significant fraction of TADs identified through wavelet analysis cannot be confirmed as GWs due to various factors, such as excessive noise leading to overly scattered Monte-Carlo results.To facilitate automation and XU ET AL.  enable statistical analysis, we have established criteria in the code for assessing whether an observed TAD event qualifies as a GW.In Table 1, we present the criteria that must be satisfied for an event to be classified as a GW.This approach aids in streamlining the process and facilitating comprehensive statistical investigations.The flowchart in Appendix E shows the data processing procedure.

Results and Discussion
We now extract TADs and identify GWs in each orbit of GOCE (November 2009 to October 2013) and in CHAMP (2004CHAMP ( -2007) ) during quiet time using the process described in Section 3. By "quiet time" we mean GW events that occurred when both the concurrent Kp and the average Kp over the past 6 hr were less than 3 in order to minimize the waves from geomagnetic disturbances.Figure 5a shows the number of TAD (gray) and GW (green) events we identified from the GOCE data during each month in 2011.During this year, the total number of GOCE TAD and GW events are 206,225 and 34,378 respectively; therefore, only ∼17% of the TAD events were GW events with relatively high confidence.For the GOCE observations during 2009-2010 and 2012-2013, the ratio of GW events identified from TAD events are similar (∼17%, not shown).Although some of the TADs are not GWs (e.g., terminator waves, etc.), other adverse conditions that prevent us from identifying GWs include but are not limited to: (a) The wavelet spectrum of a TAD wave packet can be very wide, thus making it difficult to determine the correct values of the dominant λ track and amplitudes for that wave packet.(b) The noise level from the red-noise model may have been overestimated for some GW events, yielding GW solution below the confidence threshold.(c) When the solution of a GW is close to ϕ ∼0° or 180°, it is relatively easy to measure waves in the density but almost impossible to extract cross-track wind perturbations because the amplitude is quite small.(d) When a solution has ϕ ∼ 90° or 270°, it is challenging to derive reliable GW parameters (see last row of Table 1).Therefore, the low percentage of GWs extracted from the TADs shows the limitation of our method, and should be interpreted as "at least 17% of TAD events observed in GOCE are GWs." Figure 5b shows the analogous results for CHAMP data for the year 2006.For this year (and other unillustrated years), approximately 11% of the observed TAD events are identified as GW events.This lower GW occurrence in CHAMP is potentially attributed to two factors: (a) CHAMP's observations have more noise than GOCE due to less precise accelerometers and lower density (higher altitude), thereby contributing to larger uncertainties in solved parameters following the Monte-Carlo tests; (b) the lower percentage of GWs may be geophysics due to the higher For some TAD events that are not GWs, it is impossible to have a GW solution because there is no PS or AR values in "lookup tables" that are close to PS or AR values derived from observation, but our program still returns a solution.We set these criteria to prevent those cases from being counted as solvable GWs.Unperturbed solution ∉ (70°, 110°) and ∉ (250°, 290°) The contours of -G arg and G mod at ϕ near around 90° or 270° is very dense (Figure 4), so it is easy to get a ϕ solution with large uncertainty.Additionally, according to equation λ H = λ ψ |cosϕ| (Figure 3), even a small uncertainty in ϕ can induce a large uncertainty in λ H when ϕ is around 90° or 270°.Figure 6 provides a comprehensive overview of the GW characteristics during quiet times as observed by GOCE (left column) and CHAMP (right column).In scatter plots in the top two rows in Figure 6, the c IH −τ Ir −λ H relationship generally conforms to GW dissipative theory (Vadas & Crowley, 2010).The first row shows c IH versus λ H . Two populations of GWs are apparent: the first population (with the largest # of GWs) has c IH > 300 m/s and λ H ∼ 200 km to thousands of kilometers while the second population has c IH < 300 m/s and λ H < 500 km.Thus the first population contains medium (100 km <λ H < 1,000 km) to large-scale (λ H > 1,000 km) GWs.Note that for each population, c IH is proportional to λ H , similar to previous results (Vadas & Crowley, 2010).The second row shows a scatter plot of τ Ir versus λ H .The two populations are again apparent here, although both populations have τ Ir ∼ 10−50 min.The third to fifth rows show the number of GWs as a function of λ H , c IH , and τ Ir respectively, with histograms in red and Gaussian fits in blue.There are many more GW events in the first than in the second population.The first population has a peak horizontal wavelength, intrinsic horizontal phase speed and intrinsic period of λ H = 585 km, c IH = 688 m/s and τ Ir = 14.3 min for the GOCE GWs and λ H = 655 km, c IH = 768 m/s and τ Ir = 14.4 min for the CHAMP GWs, respectively.The second population has a broad distribution of λ H = 100-400 km and peaks at c IH = 100-250 m/s with smaller (but still significant) amplitudes as compared to the first population in the GOCE data.From the CHAMP data, however, the amplitude of the second population is quite small as compared with that of the first population.This implies that the GWs in the second population have dissipated between 250 and 350 km; this is the expected behavior, because the GWs in the first population (having c IH < 300 m/s and λ H < 500 km) dissipate rapidly in the vicinity of (or below) z ∼ 250 km In the third row of Figure 6, the distribution of the horizontal wavelengths λ H of the GWs in the first population approximate a Gaussian distribution, with most of the GWs having λ H ∼ 400-800 km.R. F. Garcia et al. (2016) noted that most of the GOCE GW signal occurred in a spectral range above 8 mHz, corresponding to a maximum horizontal wavelength along the satellite track of around 1,000 km.Since their "horizontal wavelength" in that statement is equivalent to λ track in our study, because λ H ≈ λ track * |cosϕ| (where we have simplified by ignoring the Doppler's effect such that λ track ≈ λ ψ ), and since |cosϕ| ∼ 0.71 for a typical TAD in our study (since our typical TAD has ϕ ∼ 45 °), their statement is equivalent to λ H ≈ <710 km.This agrees with our results shown in Figures 6c and 6h.

General distribution
The horizontal phase speeds of the GWs in the first population are c IH ∼600-800 m/s (GOCE) and 650-850 m/s (CHAMP).In a study by Bruinsma and Forbes (2009), an investigation was conducted on 21 TAD events originating from the auroral zone and propagating across multiple CHAMP orbits into the opposite hemisphere.Their findings indicate an average speed of 646 ± 122 m/s for all TADs, which agrees with the phase speeds of most of the GWs resolved in our analysis.The horizontal phase speeds of these GWs are so large that they could not have propagated from the lower/middle atmosphere.This is because a GW can only propagate in the lower/middle atmosphere if c IH ≤ 0.9 c s , where c s is the sound speed (c s ≃ 310 m/s) in the lower/middle atmosphere (Vadas et al., 2019).Therefore, no GW with c IH ≥ 300 m/s could have propagated from the lower/middle atmosphere to GOCE/CHAMP altitudes, while the GWs in the second population could have propagated from the lower/middle atmosphere to GOCE/CHAMP altitudes.Therefore, we speculate that the GWs in the second population could be mainly primary GWs from high-frequency lower/middle atmospheric sources (such as deep convection), while the GWs in the first population are likely 1) higher-order GWs from the dissipation of GWs from the lower/middle atmosphere (e.g., from deep convection, the polar vortex, orographic forcing, jet and fronts (Becker, Goncharenko, et al., 2022;Becker & Vadas, 2020;Vadas & Becker, 2019)) and/or 2) GWs from Joule heating in the thermosphere (Hocke & Schlegel, 1996).Secondary GWs from Tropical Storm Noel, determined via reverse ray-tracing, were observed at z∼280 km at Wallops Island, VA as TIDs on 30 October 2007 with λ H = 100-2,000 km, c IH = 140-600 m/s and τ Ir = 20-60 min (Vadas & Crowley, 2010).These GWs are consistent with our first and second population results.Waldock and Jones (1986) observed medium-scale traveling ionospheric disturbances (MSTIDs) over the UK using HF Doppler sounding technique.They found that the observed MSTIDs had horizontal phase speeds of c H ∼ 50−350 m/s (with a peak at c H ∼ 140−150 m/s) and periods of τ r ∼ 10−55 min (with a peak at τ r ∼ 15−20 min).These GWs agree well with the second population of GWs observed by GOCE in Figure 6.
Using the parameters we identified for the GOCE GW events, we present the statistical results during four different seasons for 0-12 LST (Figure 7) and 12-24 LST (Figure 8).The first row of Figures 7 and 8 illustrates the global distribution of 10° × 10° lat-lon gridded mean GW intensities for various seasons.These maps show clear seasonal and LST dependence of the GW intensity.Note that the gray shading (mostly at low latitudes) display pixels containing less than 3 GW events.The summer polar regions exhibit the strongest GW intensities.In addition, at 0-12 LST, moderately strong GW intensity is observed at the summer mid-latitudes (60°S-30°S in Figure 7a and 30°N-60°N in Figure 7g).Since there is no clear geographic concentration for these GWs (i.e., no apparent concentration at specific longitudes), it is not apparent whether these mid-latitude GWs originate from  In contrast, at 12-24 LST, the GW intensities are not significantly strong at summer mid-latitudes.However, their distribution closely aligns with the location of deep convection (e.g., the African sector and North America during the local summer season in Figure 8).We explain the LST dependence of GWs from deep convection below.During the southern hemisphere wintertime for any LST, there is an increase (over the background) of GWs over and near the Southern Andes and Antarctic Peninsula at 60°S-30°S in Figures 7 and 8g (Trinh et al., 2018;Xu et al., 2021).In addition, we note that the so-called "background" distribution is quite wide in latitudinal extent (to ∼20°S) and occurs for most longitudes, thereby suggesting that the polar vortex is also a source of higher-order GWs in the winter southern thermosphere at mid to high latitudes (Becker, Vadas, et al. (2022) explored these GWs using a model and data in the northern winter thermosphere).During the northern hemisphere winter at 30°N-60°N, there is a significant increase (over the background) of GWs over the Continental United States and Greenland at 0-12 LST (Figure 7a) and even more so at 12-24 LST (Figure 8a).In addition, the "background" here also spans a wide latitudinal extent (to ∼20°N) and also covers most longitudes, thereby suggesting that the polar vortex is also a source of higher-order GWs in the winter northern thermosphere at mid to high latitudes (Becker, Vadas, et al., 2022;Vadas, Becker, Bossert, et al., 2023).Thus our observations suggest that (a) orographic forcing is an important source of higher-order GWs in the winter hemisphere, and (b) the polar vortex might be an important source of higher-order GWs in the winter hemisphere.Note that although Becker and Vadas (2020) modeled the GW hotspots during the northern and summer winter thermosphere, they attributed those hotspots as due mainly to tertiary GWs from orographic forcing; however, it is possible that some of those GWs may have instead been higher-order GWs from the polar vortex.
The second row of Figures 7 and 8 depicts the statistical analysis of GW propagation direction ξ, and the third row shows relative density amplitude (  Amp. of ρ ).The GWs in the winter and summer polar regions clearly have the largest amplitudes.We now summarize important similarities between Figures 7 and 8: 1. GW hotspot patterns occur near both geographic poles with boundaries parallel to the geomagnetic latitudes 60°S and 60°N (purple) in the first rows, especially in the polar summer season.2. The strongest GW intensities occur in the polar summer regions (first row).Additionally, the mid-and low-latitude regions (60°S-60°N) generally exhibit stronger GW activities during solstice seasons (DJF and JJA) than in the equinox seasons (MAM and SON).The low-latitude regions (30°S-30°N) typically exhibit a very low number of TAD events that are observable/solvable with our method.3. The high latitudes have significantly higher numbers of GW events compared to mid and low latitudes (second row).4. Most of the GWs have density amplitudes within the range of 0.5%-5 % (third rows).Intriguingly, all density amplitude histograms exhibit a log-normal distribution (third row).While prior research by Bezděk ( 2007)  2012) observed a log-normal distribution consistent with the ratio of accelerometer-derived densities and empirical model densities, they attributed this pattern to the multiplication of random errors, such as the satellite frontal area and drag coefficient (Doornbos, 2012).Given the ubiquity of log-normal distributions in natural phenomena (Limpert et al., 2001), unraveling the specific underlying characteristics of GWs in the thermosphere, as observed here and in other thermospheric perturbations, may warrant further inquiry to elucidate the log-normal distributions observed in these earlier investigations.
Noticeable differences also occur between Figures 7 and 8, especially in the low-and mid-latitude regions (60°S-60°N): 1.In the tropical regions, the GWs at 12-24 LST are stronger than those at 0-12 LST (first rows of Figures 7  and 8), although neither are as intense as at higher latitudes.This LST dependence likely indicates that these GWs are primary or higher-order GWs from deep convection, since deep convection has a strong local time dependence (strongest during the afternoon through local time midnight) because direct solar heating (and evaporation of water) is generally required for the atmosphere to become convectively unstable (Cotton & Anthes, 1992).Our results are consistent with Bruinsma and Forbes (2008), who found a strong LST dependence of the quiet time TADs identified from CHAMP data, with the largest amplitudes occurring for λ track ∼1,400-2,400 km at 16:00-4:00 LST.The tropical GWs at 12-24 LST are mainly concentrated near the Pacific-America and Africa sectors, corresponding to two of the peaks in the "three-peak longitudinal structure in tropics" reported in H. Liu et al. (2017).(These 3 peaks are likely related to the 3-peaked occurrence and strength of deep convective plumes in tropical regions (Vadas et al., 2014).)However, we did not observe the weakest peak from Liu's paper, which is located in the Asian-Maritime Continent sector.The absence of these GWs in our results can be attributed to their small amplitudes in both density and cross-track wind, falling below our required 95% confidence threshold using the red-noise model.This discrepancy suggests that our assumed noise model used for GW identification may be overly conservative for identifying GWs. 2. In the second row of Figure 7, the blue curve (0° < ξ < 90°) is generally the highest, and the red curve (90° < ξ < 180°) is the second-highest component for 0-12 LST.This indicates that most GOCE GWs have a propagation direction with an eastward component during this period.However, in the second row of Figure 8, either the orange curve (180° < ξ < 270°) or the purple curve (270° < ξ < 360°) is highest, suggesting that most CHAMP GWs have a propagation direction with a westward component during 12-24 LST. 3. Figures 7a and 7g demonstrate that at mid-latitudes, the summer hemisphere exhibits stronger GWs than in the winter hemisphere during 0-12 LST.Conversely, Figures 8a and 8g illustrate that at mid-latitudes, the winter hemisphere exhibits stronger GWs than in the summer hemisphere during 12-24 LST.
Further discussion on the variation of GW propagation directions and interhemispheric differences will be provided in the text related to Figures 11 and 12.
Figures 9 and 10 show the same results as Figures 7 and 8, except for the CHAMP GWs during January 2004 to December 2007 (z∼370 km).Distinct similarities and noteworthy differences between the CHAMP and GOCE GWs (Figures 7 and 8) occur.Similarities between the GOCE and CHAMP results are summarized as follows: 1. Across the same LST range and seasons, the spatial distribution of the GOCE and CHAMP GWs are similar.
As for the GOCE GWs, the polar regions consistently demonstrates higher GW intensities and event counts for the CHAMP GWs than at mid and low latitudes (first and second rows of Figures 9 and 10).2. Similar to the GOCE GWs, the CHAMP GWs also exhibit predominant eastward and westward propagation directions during 0-12 LST and 12-24 LST, respectively (second rows of Figures 9 and 10).3. Strong GW hotspots with boundaries roughly aligned with the geomagnetic latitudes of 60°S and 60°N (purple) are observed during the polar summer seasons in both GOCE and CHAMP results.
The third point suggests that GW events with substantial amplitudes are attributable to Joule heating, even during geomagnetic quiet time.This appears to contrast with Xu et al. (2021), which found that quiet time thermospheric TADs observed by GOCE and CHAMP during austral winter were predominantly induced by orographic waves and related higher-order GWs, rather than by geomagnetic activity.However, the definition of "TADs" in Xu et al. (2021) encompassed along-track perturbations spanning 160 km < λ track < 2,100 km, regardless of their density amplitude.In contrast, our current study employs a filtering mechanism that may exclude GWs with relatively smaller amplitudes that do not surpass the 95% threshold defined by the red-noise model.This filtering approach could potentially skew our results in favor of larger-amplitude GWs from geomagnetic activity over the polar regions, thereby de-emphasizing those somewhat smaller-amplitude GWs over the Southern Andes region.Additionally, while Xu et al. ( 2021) considered TADs with λ track ranging from 160 to 2,100 km, our study includes GWs with λ H from 50 to 3,000 km.It is possible that the quiet time aurora generates sizable large-scale GWs (e.g., λ H > 1,500 km), which might have been filtered out in Xu et al. (2021) but appear as hotspots in our findings.Therefore, the results of our study do not contradict the conclusions from Xu et al. (2021).As non-dissipating GWs propagate upward, their amplitudes amplify exponentially in altitude (Hines, 1960).The difference in the second point might stem from the dissipation of tropical convectively excited GWs and the generation of higher-order GWs at z∼180-200 km, which have more exponential amplitude growth when reaching CHAMP altitudes (Vadas & Liu, 2009, 2013).Or it may be due to the interhemispheric propagation of larger-scale polar-sourced GWs (Bruinsma & Forbes, 2009).For the difference in the third point, further exploration of multiple factors is provided in Figure 12 and the accompanying text.

Differences between GOCE and CHAMP results shown in
We now investigate the propagation direction of the GOCE and CHAMP GWs as a function of LST.In particular, significant differences are seen in the GW propagation direction ξ between 0 and 12 LST and 12-24 LST in the middle rows of Figures 7-10; for both GOCE and CHAMP, most GWs have propagation directions with an eastward and westward component during 0-12 LST and 12-24 LST, respectively.In Figure 11 we show a more complete relationship between ξ and LST. Figure 11a shows a scatter plot of ξ and LST for all GW events in the northern hemisphere at 0°-90°N.To see the dominant GW propagation direction more clearly, we divide ξ into  four quadrants (blue, red, orange, purple for ξ ∈ (0°,90°), (90°,180°), (180°,270°), (270°,360°), respectively) and bin hourly in 24 LST bins. Figure 11b shows the percentage of GW events having azimuths within these 4 different ξ quadrants as a function of LST.Figures 11c and 11d shows ξ versus LST for GWs observed by GOCE in the southern hemisphere at 90°S-0°.As mentioned in Section 2, the LST of GOCE is 6:13-7:33 during the descending nodes and 18:13-19:33 during the ascending nodes, so the LST of GOCE's measurements are highly dependent on latitude.Some LST dependence of ξ is visible for certain LSTs.For example, around 6 LST in the northern hemisphere, the percentage of GWs with 0° < ξ < 90° (blue) decreases, while that for GWs with 90° < ξ < 180° increases (see Figure 11b).The trend is opposite right after 6 LST in southern hemisphere (see Figure 11d).However, it is difficult to determine how the GOCE GW ξ depends on LST due to the limited LST range.
Figures 11e−11h show the corresponding results for the CHAMP GWs.Because CHAMP orbit drifted rapidly, CHAMP covered all LST.However, because of CHAMP's nearly 90° inclined orbit (87°), there are subtle yet discernible gaps in the data in Figures 11e and 11g near ξ ∼ 0°, 90°, 180°, and 270°.These characteristic gaps closely align with the gaps from our analysis method, since our method does not solve for GWs that propagate parallel or perpendicular to CHAMP's mostly north/southward tracks.From Figure 11f, during LST of 1-7 hr/7-13 hr/13-20 hr/20-1 hr, the GW azimuths are 0° < ξ < 90° (blue)/90° < ξ < 180° (red)/180° < ξ < 270° (orange)/270° < ξ < 360° (purple), respectively.In contrast, Figure 11h shows that during LST of 2-7 hr/7-13 hr/13-19 hr/19-2 hr, the GW azimuths are 90° < ξ < 180° (red)/0° < ξ < 90° (blue)/270° < ξ < 360° (purple)/180° < ξ < 270° (orange), respectively.Therefore, the GW propagation direction in northern and southern hemispheres clearly show a strong clockwise and counterclockwise rotation in a local day, respectively.This is the expected result, because the GWs that can best survive dissipation and critical level filtering propagate in a direction opposite to the background wind (Becker, Goncharenko, et al., 2022;Becker, Vadas, et al., 2022;Fritts & Vadas, 2008;Waldock & Jones, 1984).Since the background thermosphere wind (due mainly to the diurnal migrating tide from extreme ultraviolet EUV heating at these altitudes) rotates clockwise (counterclockwise) in the northern (southern) hemisphere (Cowling et al., 1971), clockwise (counterclockwise) rotation of the GWs is expected to occur in the northern (southern) hemisphere.Indeed, clockwise rotation of the GW propagation direction was observed using HF Doppler of MSTIDs over the UK (Waldock & Jones, 1986).However, this work found the estimated angle between the TID azimuth ξ to be ∼130°-140° clockwise with respect to the wind direction (the propagation direction of TIDs lag the anti-wind direction by 30°-40°).In addition, TIDDBIT observations over Wallops Island, VA, USA also found a clockwise rotation in time with a relative angle between ξ and the background wind direction of ∼120°-160° (based on our interpretation of the data shown in Figures 9a and 10a of Crowley and Rodrigues (2012), although they state the relative angle to be ∼ 90° in the text).Crowley et al. (1987) observed a counterclockwise rotation of wave azimuth in the southern hemisphere.Importantly, Figures 11f and 11h show that the times when the GWs propagate mostly southward occurs at ∼13:30 LST.This is the approximate time (or ∼1 hr later) than when the wind from the migrating diurnal tide is northward (Becker, Vadas, et al., 2022).Thus we find a relative angle between ξ and the background wind direction of ∼165°-180°.Know that the azimuth of a TID generated by a GW via ion-neutral collisions is the same as the azimuth of this GW for a single ion species (z > 180 km) (Nicolls et al., 2014).Therefore, our results agree very well with these previous observations, thus providing a good validation of our method.
Previous studies have shown that GW activity in the lower atmosphere significantly influences the thermosphere/ionosphere during quiet times, including deep convection during the summer (H.Liu et al., 2017;Trinh et al., 2018;Vadas & Liu, 2013;Vadas et al., 2014), and higher-order GWs from wind flow over orography (Becker & Vadas, 2020;Vadas & Becker, 2019) and from the polar vortex in the winter (Becker, Vadas, et al., 2022;Frissell et al., 2016).However, a quantitative understanding of the relative importance of the sources of GWs "from below" is still lacking.To address this, we investigate the monthly changes of the GW amplitudes.In Figure 12, we show the monthly multiplicative mean (μ*) of the relative density amplitude (  Amp.of ρ ) of GOCE and CHAMP GWs for specific latitude ranges.The first striking observations is that the GW amplitudes at any LST and latitude range show a clear semi-annual variation for both GOCE and CHAMP, being maximum (minimum) at the solstices (equinoxes) in both the northern and southern hemispheres.However, this semi-annual variation is much more important for the GOCE GWs than for the CHAMP GWs, with the solstice maxima being of similar amplitudes.In the CHAMP data, however, the annual variation is much more important, with the JJA maxima being much larger than the DJF maxima for both the northern and southern hemisphere.This difference is presumably due to the dissipation of GWs from molecular viscosity between the different satellite altitudes XU ET AL.

10.1029/2023JA032078
21 of 35 of 270 and 370 km, and would occur if the DJF GWs have slower horizontal phase speeds than the JJA GWs (Vadas, 2007).
During 0-12 LST, the GOCE GW intensity is consistently higher at 60°S-60°N (and 0°-60°N) during JJA than at 60°S-60°N (and 60°S-0°) during DJF, with the largest contribution being from 0° to 60°N (blue curve) (see Figure 12a).Thus GWs during 0-12 LST at low-and mid-latitudes are strongest during JJA, which coincides with the peak of the precipitation at 0°-20°N (Cullens et al., 2022).It also coincides with the time of the largest GW peak over the Southern Andes in GW observations near mid-stratospheric (L.Hoffmann et al., 2013) and stratopause (Xu et al., 2023).During 12-24 LST, however, the GWs generally exhibit larger amplitudes during DJF, with the largest contribution being from 0° to 60°N (blue curve) (see Figure 12b).This indicates that the wintertime GWs during 12-24 LST have a stronger impact on the thermosphere than those in the summer thermosphere.For comparison, we also plot the μ* values of GWs for each hemisphere and for the entire globe.As shown in the lower two panels, when GWs in high-latitudes (90°S-60°S and 60°N-90°N) are included, the summer hemisphere consistently exhibits stronger GWs, both for GWs in 0-12 LST (Figure 12c) and 12-24 LST (Figure 12d).This means that the amplitudes of GW events in the polar summer region are generally stronger.The overwhelming numbers of GW occurrence in summer polar region effectively shift μ* values of the whole summer hemisphere higher.Figures 12e and 12f show the monthly variation of the CHAMP GW density amplitudes.The GWs at any LST again are strongest during JJA, which coincides with the peak of the rainfall at 0°-20°N.Cullens et al. (2022) found that the amplitudes of the GWs observed by ICON-MIGHTI at 0°-40°N exhibited semi-annual variation in the mesosphere-lower-thermosphere (MLT), with maximum (minima) at the solstices (equinoxes).That study also showed that the GW zonal wind perturbations, u′, at 200 km altitude exhibited mainly an annual variation with a weaker semi-annual variation, with the major (minor) maxima occurring at JJA (DJF) and with the minima occurring at the equinoxes (their Figure 2a).The difference between their results and our GOCE GWs (shown in Figure 12) is that their DJF MIGHTI u′ peak is 40% weaker than their JJA peak.This difference likely occurs for the following reason.The dominant sources of the thermospheric GWs during DJF at 0°-40°N (as observed by ICON-MIGHTI) are likely the polar vortex and orographic forcing.Because these GWs would have had to propagate large distances meridionally in order to be observed by MIGHTI at 0°-40°N, most would have had larger meridional than zonal components of their propagation directions, which is equivalent to having larger meridional (than zonal) velocity perturbations.Because Figure 2a in Cullens et al. (2022) only showed u′, that figure would have underestimated the total GW amplitude (i.e.,  √ ( ′ ) 2 + ( ′ ) 2 ) during DJF by  √ 2 ∼1.4,which is a false DJF peak reduction of >30%.Taking this factor into account, our GOCE GW results agree well with the ICON-MIGHTI results.
Overall, the results in GOCE and CHAMP are similar to the results observed by HF Doppler at a similar altitude (Figure 6a of Waldock and Jones (1986), which shows the GW frequency peaks in February and July).We note that the GW amplitudes measured at mid to high latitudes using a meteor radar in the MLT also exhibit a semi-annual variation, with maxima at the solstices (e.g., Figures 11 and 12 of P. Hoffmann et al. (2010)).That behavior is due to the annual variation of the mean wind in the MLT region over mid-and high-latitude regions, which goes through zero (and reverses sign) at the equinoxes (e.g., Figure 4b in R. R. Garcia et al., 1997;P. Hoffmann et al., 2010).When the mean wind is close to zero, the GWs propagating in all directions break and deposit their momentum and energy in the MLT.However, during solstice, those GWs propagating opposite to the mean wind propagate much higher before breaking or dissipating in the thermosphere, thus projecting their semi-annual variation onto the higher-order GWs.Thus, the semi-annual variation observed in Figure 12 here strongly suggest that these GOCE and CHAMP GWs are mainly higher-order GWs from the dissipation of primary GWs from GW sources in the lower/middle atmosphere.

Conclusions
In this study, we investigated the thermospheric GWs observed by the GOCE and CHAMP satellites.GWs are essential for providing variability to the thermosphere, generating TIDs, redistributing energy and momentum throughout the Earth's atmosphere, and influencing the dynamics of the thermosphere over a range of spatial and temporal scales.Recent research has shown a strong link between lower atmospheric dynamics, such as polar vortices, and the characteristics of GWs in the thermospheric/ionospheric region.XU ET AL.

10.1029/2023JA032078
22 of 35 Although various methods have been employed to study thermospheric GWs, in-situ spacecraft missions equipped with precise accelerometers offer unprecedented accuracy and resolution.Our study leveraged this advantage to determine the intrinsic characteristics of GW events observed as TADs along the satellite tracks.We used a wavelet analysis and the red-noise model to analyze each TAD event and determine if it was a thermospheric GW event.The application of wavelet in this current study was based on 2 assumptions.(a) The spectrum of noise is continuously distributed and can be modeled by the red-noise model.In this simplified model, noise is assumed to have a power (variance) spectral density that decreases as the scale decreases.(b) The GWs that superimposed on the red-noise sequence covers a narrow spectrum instead of a broad dispersive spectrum, so these GWs can be well represented by sinusoidal waves that are superimposed on the red-noise series.Our results showed that those assumptions were generally satisfied but were somewhat conservative for diagnosing TADs in general, because we saw many small-scale TADs that were filtered out when we applied the 95% confidence level as the threshold to extract TADs for our method.
In sum, we analyzed the quiet-time (Kp < 3) TADs from GOCE during November 2009 to October 2013 (at z ∼ 270 km) and from CHAMP (at z ∼ 370 km) during 2004-2007.Once we extracted the density and cross-track wind perturbations for the TADs, we used the GW dissipative polarization and dispersion relations, combined with the background thermosphere, to derive the intrinsic properties of some of the GW events (after applying the selection criteria given in Table 1), such as their propagation directions, intrinsic periods, and horizontal wavelengths.Despite the above challenges, our methodology provides valuable insights into the statistical properties of GWs.Here are the findings: 1. We found that ∼17% of TADs in GOCE and 11% of TADs in CHAMP were confidently identified as GW events.2. The peak horizontal wavelength, intrinsic horizontal phase speed and intrinsic period were λ H = 585 km, c IH = 688 m/s and τ Ir = 14.3 min for the GOCE GWs and λ H = 655 km, c IH = 768 m/s and τ Ir = 14.4 min for the CHAMP GWs, respectively.The increase of c IH with height is consistent with GW dissipation from molecular viscosity.Such fast phase speeds show that most of these GWs were generated in the thermosphere, not in the lower/middle atmosphere (for which it is necessary that those GWs have c IH <0.9*310 m/ s∼ 280 m/s).3. We also find a second population of GOCE GWs with λ H = 100-400 km, c IH = 50-250 m/s and τ Ir = 10-50 min.
These slower GWs were also present in CHAMP but much less frequently.It is possible that many of these second-populations GW were from the lower/middle atmosphere.
Next, we provided global morphology maps of the thermospheric GWs (instead of global TAD maps) from GOCE and CHAMP for the first time.The findings derived from these maps are: 1. High latitudes consistently exhibit higher GW intensities, with polar summer regions showing the strongest GW amplitudes.Although a significant amount of the polar GWs in thermosphere may be higher-order GWs from lower/middle atmospheric processes (such as orographic forcing, the polar vortex, and deep convection), some of our diagnosed GWs may also be generated by Joule heating, even though we attempted to minimize these latter GWs by restricting Kp < 3 for allowed events.2. We also found that higher-order GW hotspots occurred during the winter thermosphere at mid to high latitudes, likely from primary GWs from orographic forcing and the polar vortex.We found that the GOCE GW amplitudes exhibited semi-annual variations, with maxima in January and July, and minima at the equinoxes.In contrast, the CHAMP GW amplitudes exhibited mainly annual variation (with a weaker semi-annual variation), with the major (minor) maximum in July (January), and with minima at the equinoxes.These seasonal variations show that a significant fraction of the GOCE and CHAMP GWs (identifiable by our method) are likely not from Joule heating, but are higher-order GWs originating from the dissipation of primary GWs generated in the lower/middle atmosphere.This dissipation is shaped by seasonally varying background winds during GW propagation and primary GW sources.
Finally, we found that the average propagation direction of the CHAMP GWs exhibited a clear diurnal cycle, with clockwise (counter-clockwise) occurring for GWs in the northern (southern) hemisphere, in agreement with previous observations of TIDs.In addition, we found that the CHAMP GWs had equatorward propagation occurring at ∼13:00 LST.This suggests that the predominant propagation direction of CHAMP GWs is approximately opposite to the background wind direction.According to Torrence and Compo (1998, hereafter TC98), for a discrete signal sequence x n , the discrete Fourier transform (DFT) of ) where n = 0…N−1 is the index of x n and k = 0…N−1 is the frequency index.In practice, x n usually represents a time series (e.g., in TC98), and each increment in n corresponds to the time step δt.However in the current study, x n represents a series along the satellite track with length step δx (which is the ground-based distance that satellite travels in 10 s).The operator (^) denotes the Fourier Transform.The complex wavelet spectrum W(n, s) is defined by the convolution theorem: where ω k is the angular frequency defined by ω k = sgn(N/2 − k)(2πk/Nδx).The wavelet function    () is given by Equation 6 in TC98.We designate the Morlet wavelet defined in their Table 1 to be wavelet base function   0 .The (*) indicates the complex conjugate.Following Section 3f in TC98, the wavelet scale s is written as fractional powers of two and in the unit of the spatial resolution δx along the satellite track: where s 0 is the smallest resolvable scale.In this study, we use δj = 0.25 and s 0 = 2δx for this paper.J determines the largest scale of s j within the length of the input series x n .
One of our assumptions is that when no TAD signal is present in the GOCE or CHAMP data sets, both the variations in the density and cross-track wind along the orbit can be modeled by red noise.Following TC98, a simple model for red noise is the univariate lag−1 autoregressive [AR(1), or Markov] process: where α is the assumed lag-1 autocorrelation, x 0 = 0, and z n is Gaussian white noise.Following Gilman et al. (1963), the discrete Fourier power spectrum at frequency index k (here k = 0…N/2), after normalizing with the factor 1/σ 2 (where σ 2 is the variance of x n ), is: The meaning of Equation B5 is: noise has a power spectral density that decreases as the scale decreases.Since each    at a fixed k is normally distributed in the red noise model, then the wavelet coefficients (the band-passed inverse Fourier components) should also be normally distributed (TC98).Therefore, when the wave signal is absent or subtracted from the series x n , W(n, s) is normally distributed in the complex plane at (n, s), and the wavelet power spectrum |W(n, s)| 2 of red noise is chi-square distributed: where   2 2 is the chi-square distribution with two degrees of freedom.According to the chi-square distribution , we say that it is a real wave signal (instead of noise) in the observations with confidence higher than 95%, because the probability of this peak to be noise is less than 5%.
As shown by the in Maraun and Kurths (2004, see their Figure 2) and pointed out by Y. Liu et al. (2007), the wavelet power spectrum W(n, s) calculated by TC98's method is good when the input is dominated by stochastic (noisy) processes, but has a bias toward larger scales s when the noise input is superposed by monochromatic XU ET AL.
Following C16, we derive the unbiased wavelet spectrum  mod( ) using   mod() when we assume the input has no noise: Since the complex wavelet spectrum W(n, s) is expected to have a 2D normal distribution in the complex plane given by Equation B5, the expected uncertainty (conventionally referred to as "half-length of error bar") in amplitude can be quantified by the standard deviation of either  ℜ{mod( )} or  ℑ{mod( )} , which is For TADs with high peaks in amplitude spectrum (i.e., when  Amp  () ≫ SDAmp( ) , which brings high confidence to say it is not noise), the corresponding uncertainty in phase can be represented by Note that the counterpart of the "equivalent Fourier period" in TC98 is the along track wavelength (λ ) in our study, and thus gives a value of λ track = 1.03s for the Morlet wavelet for the nondimensional frequency ω 0 = 6.Therefore, variables that can be written as functions of (n, s) can be written as functions of (n, λ track ) as well.Based on the expectation of W mod defined by Equation B9, the distribution of W mod in the complex plane is: Appendix C: Deriving the PS and AR of Perturbations of  ρ and  û Along the Satellite Track Direction (ψ) as Functions of ϕ and λ z Following Vadas et al. (2019), the complex, dispersion relation for acoustic waves and GWs damped by molecular viscosity and thermal diffusivity is a function of wave number vectors k H = (k, l) and m and the complex intrinsic frequency ω I : )) ] = 0. (C1) XU ET AL.
Respectively.Here k = (k, l, m) = 2π/(λ x , λ y , λ z ) is the full wave vector.The background characteristics in equations ( C1)-(C3), such as N B (buoyancy frequency), c s (sound speed),   (density scale height), ν (kinematic viscosity,   = ∕ ,  is the molecular viscosity), γ (ratio of specific heats, γ = c p /c v , in which c p and c v are isobaric and isochoric specific heat respectively), etc., are determined directly or indirectly by gravity and background density using the NRLMSIS2.1 model, while α and   are functions of these characteristics and k H and m.Finally, a in Equations C1-C3 is set to be 1 to include the bulk viscosity in addition to the shear viscosity in the viscous stress tensor.Please refer to Section 4.3 in Vadas et al. (2019) for further information.
By dividing Equations C2 and C3, the vertical wind perturbation  w on the left-hand-sides (LHSs) of Equations C2 and C3 cancel out.The ratio of the right-hand-sides (RHSs) of Equations C2 and C3 is still a function of λ H and λ z .Therefore, the ratio of RHSs of Equations C2 and C3 can be expressed as functions of the horizontal and vertical wavelengths.We use F to denote this function: (C11) Although G in Equation C10 is a function of three independent variables (λ H , −λ z , ϕ), λ H is dependent on ϕ.
(The function λ H with respect to ϕ, i.e., k H = 2π/λ H = H(ϕ, −λ z ), is derived in Appendix D. See the derivation of Equation D7) Therefore, G can be reduced to a function of unknown independent variables that Equation C10 can be written as  ( ρ∕ ũtrack) = ( −) .Since the modules and argument of its LHS  ( ρ∕ ũtrack) can be extracted from amplitude and phase spectra using unbiased wavelet analysis, the two unknowns (ϕ, −λ z ) can be solved by using charts of PS ψ and AR ψ as "look up tables".The other GW intrinsic parameters, such as, λ H , τ Ir , ξ, etc., can be derived directly based on (ϕ, −λ z ) solutions.As shown in Figure D1, the distance between the phase lines of the GW perturbations along the direction of azimuth ψ at a fixed time t is λ ψ .The along-track wavelength of the GW perturbations observed by the satellite is λ track .If we assume the satellite velocity v sat is much larger than the ground based horizontal phase speed of a GW, then the observed wavelength along track λ track is very close to λ ψ (i.e., when v sat →∞, λ track →λ ψ ).This assumption may not be reasonable in extreme scenarios, especially when |cosϕ| is small so the satellite velocity along the GW propagation direction |v sat cosϕ| and the GW horizontal phase speed c H are comparable.Therefore, it is necessary to determine λ ψ as a function of λ track for the extracted measurements before solving for the GW intrinsic parameters.

Figure 1 .
Figure 1.The procedure to identify a traveling atmospheric disturbance (TAD) event from GOCE density perturbations and cross-track wind.(a) The positions of GOCE data points over one full circle track divided by "DA mode" during Coordinated Universal Time (UTC) 2010-07-05_22:10:20∼23:39:50 (red dots), subsolar point (yellow sun symbol) and day-night terminator (yellow curve).The red westward arrow denotes the satellite direction when GOCE starts this track.The open/ filled red circle denotes the first/last position of the data point in this track.(b) GOCE in-situ measured density ρ obs in log scale versus altitude within this orbit (black).The density scale height is derived from the slope of the linear-fitted line (red) of the log-scaled densities.(c) The black curve presents density (ρ) that are extrapolated to z = 270 km according to the measured thermospheric density scale height as a function of distance along the satellite track.The blue dashed line denotes the mean extrapolated density in this orbit at z = 270 km (   ).(d) Relative density perturbation in percent,   ′ ∕ (black, where   ′ =  −  ).The green curve denotes the background density for wavelengths ≥3,000 km.(e) Unbiased amplitude spectrum of the detrended (λ track < 3,000 km) relative density perturbation (   ′ ∕ ).The thick black contour encloses TADs with >95% confidence for a red-noise process.Within the >95% confidence contours, peaks labeled with green crosses and green indices denote locations of TADs in density.Gray cross-hatched regions indicate the "cone of influence," where edge effects are important and waves are not trustworthy.(f-g) The same as (d-e) but for cross track wind perturbation   ′ track .In panel (g), peaks labeled with orange crosses and orange indices denote locations of TADs in cross-track wind.(h) Comparison between the sum of reconstructed perturbations    (blue, λ track < 3,000 km, confidence>95%) and the original detrended   ′ ∕ (black, λ track < 3,000 km).The horizontal axis only spans from 2 × 10 4 to 3 × 10 4 km.(i) The same as (g) but for reconstructed perturbations    (orange, λ track < 3,000 km, confidence>95%) and detrended   ′ track (black, λ track < 3,000 km).(j) The reconstructed perturbation pair of , in %) and horizontal wind along the GW propagation (

Figure 2 .
Figure 2. (a) The phase shift between    and    for gravity waves (GWs) along their propagation direction (contour labels are in unit of degrees).(b) The amplitude ratio between    and    , that is,  |100 * ρ∕ ũ | for the GWs (contour labels are in unit of %*s/m).The background state of the thermosphere is determined by NRLMSIS2.1 for 5 July 2010 at 22:53:00 UTC, 79.43°S, 75.53°E, and z = 270 km.(a) and (b) are also visualization of −F arg (λ H ,−λ z ) and |F(λ H , −λ z )|, respectively, as defined by Equation 1.

Figure 3 .
Figure 3. Track of the satellite which is moving with azimuth ψ (thick black arrow), gravity wave (GW) propagation direction (thick magenta arrow with azimuth ξ), ϕ = ψ−ξ(clockwise from ξ to ψ) and GW lines of constant phase (green parallel lines) in the horizontal plane.λ H is the GW horizontal wavelength, defined as the distance between two consecutive GW phase fronts at a fixed moment.λ ψ is the distance between two consecutive GW phase fronts measured by the satellite along the track at a fixed moment (defined as λ H = λ ψ * |cosϕ|).(Note that λ ψ and λ track are different because the former quantity is the actual distance between phase fronts after taking into account the Doppler effect, while the latter quantity is the distance between peaks measured by the satellite.Details are contained in Appendix D and Figure D1.) "Sat.Q-I" through "Sat.Q-IV" denotes the 4 Satellite quadrants relative to the satellite flight direction.

Figure 4 .
Figure 4.An example of solving a gravity wave (GW) event.(a) Phase shift (unit in degrees) between    and   track for a GW detected along the satellite track with λ track = 1,087.2km as a function of ϕ = ψ−ξ.(b) The same as (a) but for the GW amplitude ratio  |100 * ρ∕ ũ | (unit in %*s/m).(a) and (b) are also visualizations of -G arg (λ H , −λ z ) and |G(λ H , −λ z )|, respectively, as defined by Equation 4. Green contours represent  PS( ρ ρtrack) = 154.3°,and orange contours represents  AR( ρ ρtrack) = 0.105, without considering any noise in the measurements.Magenta plus signs represent 1,000 Monte-Carlo solutions of    and   track that satisfy the red-noise model.The intersection of most Monte-Carlo solutions in Sat Q-II (ϕ∈(0 °, 90 °)) are considered the "unperturbed solution" (also optimal solution in this case) and are denoted by blue diamonds.(c) Occurrence frequency distribution of intrinsic GW parameters (λ H , τ Ir , ϕ, λ z ) derived from the 1,000 Monte-Carlo solutions (purple).The blue curves are Gaussian (ϕ) or log-normal (λ H , τ Ir , λ z ) fits to each histograph.(d) Data points along the satellite track (red dots), location of the observed GW (green cross), the direction of satellite motion (black arrow), and the direction of the GW propagation obtained by the peak of  in (c) (lower left).
in solution/observed PS) ∈ (1 ± 30%) AND Unperturbed AR in solution ∈ (observed AR ± 30°) (370 km).The quantification of the roles these two factors play in GW percentages within TADs necessitates further investigations targeting the non-GW TADs.

Figure 5 .
Figure 5. (a) Number of GOCE traveling atmospheric disturbance (TAD) events identified each month in 2011 (gray) alongside gravity wave (GW) events (green).The total count of TAD and GW events are 206,225 and 34,378, respectively.(b) The same as (a) but for CHAllenging Minisatellite Payload TADs and GWs during 2006, with a total of 202,213 TAD events (gray) and 22,270 GW events (green).

Figure 6 .
Figure 6.Attributes of all gravity wave (GW) events observed by GOCE (left column) and CHAllenging Minisatellite Payload (CHAMP) (right column).The total number of identified GW events is given at the top of each column.(a) Scatter plot of intrinsic phase speed c IH versus horizontal wavelength λ H ; (b) Scatter plot of intrinsic period τ Ir versus λ H ; (c), (d), and (e) represent the number of GW events as a function of λ H , c IH and τ Ir , respectively (red histograms).Gaussian fits are provided for (c-(d) (blue lines) with the μ (mean) and σ (standard deviation) given in each panel.(f-g) The same as left column (a-e) but for GWs observed in CHAMP.

Figure 7 .
Figure 7. Statistical results based on gravity waves (GWs) observed by GOCE during 0-12 local solar time during November 2009 to October 2013.Each column shows the GW events in four different seasons: DJF (December-January-February), MAM (March-April-May), JJA (June-July-August), and SON (September-October-November), respectively.(a) Global distribution of gridded mean GW intensity or  ∑  (Amp. of ρ) 2 ∕ , where  (Amp. of ρ) represents the relative density amplitude (unit %) of an identified GW event and N is the total number of identified GWs in a 10° × 10° map grid.Gray grids indicate regions with less than 3 GW events.The purple curves in maps depict the geomagnetic latitude of 60°, 0°, and −60°.(b) Number of GW events in four different propagation direction quadrants as a function of latitude.The blue curve corresponds to 0° < ξ < 90°, the red to 90° < ξ < 180°, the orange to 180° < ξ < 270°, and the purple to 270° < ξ < 360°.(c) The number of GW events as a function of the relative density amplitude  (Amp. of ρ) (red histograms).Note the logarithmic scale on the horizontal axes.Gaussian fits (blue), and the fit parameters μ* (multiplicative mean) and σ* (multiplicative standard deviation) are displayed in each panel (blue).(d-f) show the same analysis as (a-c) but for the MAM season.(g-i) show the same analysis as (a-c) but for the JJA season.(j-l) show the same analysis as (a-c) but for the SON season.
-latitude deep convection or are generated from the polar regions and propagate to mid latitudes.

Figure 8 .
Figure 8.The same as Figure 7 but for gravity waves observed by GOCE during 12-24 local solar time.
Figures 7-10 are: 1.The intensity (or relative density amplitude) of the GWs observed by CHAMP are generally larger than that observed by GOCE (see the color bars and histograms in the third rows).2. While GOCE only identifies pronounced tropical GW clusters near the Pacific-American and African sectors, CHAMP records elevated intensities and occurrences of equatorial GWs across both 0-12 LST and 12-24 LST.Specifically, the number of GW events around the equator outnumber those at latitudes of 30°S and 30°N (compare the second rows of Figures 7-10).

Figure 9 .
Figure 9.The same as Figure 7 but for gravity waves observed by CHAllenging Minisatellite Payload during 2004-2007 (0-12 local solar time).Note that the color bar range is larger than that for GOCE in Figures 7 and 8.

Figure 10 .
Figure 10.The same as Figure 8 but for gravity waves observed by CHAllenging Minisatellite Payload during 2004-2007 (12-24 local solar time).

Figure 11 .
Figure 11.(a) Scatter plot depicting the azimuthal gravity wave (GW) propagation direction ξ as a function of local solar time (LST) for GOCE GWs (2009-2013) within the northern hemisphere.Each data point represents an individual GW event.(b) The distribution of GW events across the four ξ quadrants binned as a function of LST.The blue curve corresponds to 0° < ξ < 90°, red to 90° < ξ < 180°, orange to 180° < ξ < 270°, and purple to 270° < ξ < 360°.Panels (c) and (d) adopt the configuration of (a) and (b) but for GOCE GWs within the southern hemisphere.(e-h): The same as (a-d) but for the CHAMP GWs (2004-2007).Because of the quasi-fixed local time of GOCE's orbit during 2009-2013, the LST of GOCE's measurements along orbit are highly dependent on latitude.To show this dependency, gray approximated latitude markers are added above the upper y-axis in panels (a) and (c).

Figure 12 .
Figure 12. (a-b) Time series of monthly multiplicative mean (i.e., log-normal fitted peak values μ*). of the relative density amplitude of gravity waves (GWs) (in %) at latitudes of 60°S-60°N (black), 60°S-0° (red) and 0°-60°N (blue).(a)/(b) shows results for GWs at 0-12/12-24 local solar time, respectively.(c-d) Same as (a-b) butfor GOCE GWs at all latitudes (black), in the southern hemisphere (red) and in the northern hemisphere (blue).(e-f) Same as (c-d) but for CHAllenging Minisatellite Payload GWs at all latitudes (black), in the southern hemisphere (red) and in the northern hemisphere (blue).Thick arrows at January and July indicate which hemisphere has the higher mean μ*, with red/blue arrows for the southern/northern hemisphere.

Appendix B :
Deriving the Unbiased Wavelet Spectrum W mod (n, s) and Determine the Corresponding Uncertainties.
reduce this bias, Chen (2016, hereafter C16) replaced the wavelet function    () in TC98 by a modified wavelet function   mod()  .The conversion between them is ψmod GW complex polarization relations for the density and horizontal wind perturbations (substituting    → Ã ,  ŵ → w ,    → Ã into Equations 11 and 13 of Vadas et al. (2019), since the operator ^ has already been used to denote the Fourier Transform in this paper) along GW propagation direction ξ can be expressed as two equations: Figure 2 visualizes  PS( ρ ρ ) and  AR( ρ ρ ) at a specific time (UTC 2010-07-05 22:53:00) and specific location (79.43°S, 75.53°E, z = 270 km) during the quiet-time austral winter (Kp index = 1−).To derive the PS and AR between two perturbation variables of a GW along the satellite track azimuth ψ based on (C5), the value of ϕ should be taken into consideration.When ϕ∈(−0.5π,0.5π), the PS of    and    along the satellite direction ψ is equal to  PS( ρ ρ ) , otherwise they have the opposite sign.Thus we have PS ( ρ ρ ) = sgn(cos ) * PS( ρ ρ ).(C6) Because   track = −  * sin  (as shown in Figure 3), we get Appendix D: Derivation of λ H of a GW as a Function of ϕ, Background Wind, and Satellite Movement Parameters.

Figure D1 .
Figure D1.(a) Illustration of sampling effect when  is in the range of (−0.5π, 0.5π), corresponding to scenario when the satellite overtakes the GWs; (b) The same as (a) but when  is in the range of (0.5π, 1.5π), corresponding to scenario when the satellite meets the gravity waves.(c) Diagram of function  PS( ρ ρtrack) = −arg( −) without considering sampling effect by assuming v sat = infinite for GOCE satellite; (d) The same as (c) but considering sampling effect with v sat = 7,765 m/s for GOCE satellite.

Figure F1 .
Figure F1.Methodology application using the GOCE gravity wave (GW) cases generated by Tohoku-Oki earthquake at 05:46 UTC on 11 March 2011.(a) The reconstructed perturbation pair of  10 * ρ (blue) and    (red) at λ track = 384.4km measured in track 1.The "PS" means "phase shift" in units of degrees, "AR" means "amplitude ratio" in units of %*s/m.(b) The same as (a) but for λ track = 420.7 km measured in track 2. (c) The same as (a) but for λ track = 841.5 km measured in track 3. (d) The map with GOCE data points (red dots), the location of earthquake (yellow cross) and locations of peak amplitudes measured by GOCE (green crosses).The magenta arrows denote the propagation directions and λ H of the solved GWs.Solved intrinsic parameters and uncertainties of the GWs are shown in the magenta boxes.
shows the unbiased wavelet spectrum |W mod | in °, 270 °)) in ϕ−λ z space, meaning that there are two analytical solutions if we assume no uncertainties caused by noise.Then using the red-noise model described in Section 3.1 and Appendix B, we generate 1,000 Monte Carlo values that satisfy the distribu- Percent of GW solutions in Monte-Carlo tests >60% If there are too many (>40%) Monte-Carlo tests that do not return GW solutions, then the Monte-Carlo tests are considered invalid General distribution Percent of GW solutions in one Sat.Q > 60% If GW solutions provided by Monte-Carlo tests are overly spaced-out in different satellite quadrants, then the Monte-Carlo tests are considered invalid Unperturbed solution of λ H > 50 km Lower limit for reliable λ H solutions, since the 10 s sampling cadence corresponds to ∼80 km along-track distance for GOCE and CHAMP.Ir Unperturbed solution of τ Ir > 5 min Lower limit since τ Ir > the buoyancy period τ B , which is > 5 min in the thermosphere ϕStandard deviation of ϕ < 30°ϕ

Table 1
Criterion Required to Determine if a Traveling Atmospheric Disturbance Event Is a Gravity Wave Event XU ET AL.
Zonal and meridional components of the background wind.U H : Background horizontal wind along the direction of GW propagation.azi U : azimuth of horizontal direction of U.
arg($) or $ argArgument of a complex number, array, or function $.}Real and imaginary part of a complex number, array, or function $.

table ,
− 1 ), the factor of two in the right-hand side of the equation is needed to derive the correct amplitudes.Therefore, if uncertainties are included, the corresponding expectation of the amplitude and phase at ( ) is defined as: