Satellite observations of gravity wave momentum flux in the mesosphere and lower thermosphere (MLT): feasibility and requirements

In the recent decade it became evident that we need to revise our picture of how gravity waves (GWs) reach the mesosphere and lower thermosphere (MLT). This has consequences for our understanding not just of the properties of the GWs themselves, but in particular of the global circulation in the MLT. Information on spectral distribution, direction, and zonal mean GW momentum flux is required to test the theoretical and modeling findings. In this study, we propose a constellation of two CubeSats for observing mesoscale GWs in the MLT region by means of temperature limb sounding in order to derive such constraints. Each CubeSat deploys a highly miniaturized spatial heterodyne interferometer (SHI) for the measurement of global oxygen atmospheric band emissions. From these emissions, the 3-D temperature structure can be inferred. We propose obtaining four independent observation tracks by splitting the interferograms in the center and thus gaining two observation tracks for each satellite. We present a feasibility study of this concept based on self-consistent, high-resolution global model data. This yields a full chain of end-to-end (E2E) simulations incorporating (1) orbit simulation, (2) airglow forward modeling, (3) tomographic temperature retrieval, (4) 3-D wave analysis, and (5) GW momentum flux (GWMF) calculation. The simulation performance is evaluated by comparing the retrieved zonal mean GWMF with that computed directly from the model wind data. A major question to be considered in our assessment is the minimum number of tracks required for the derivation of 3-D GW parameters. The main result from our simulations is that the GW polarization relations are still valid in the MLT region and can thus be employed for inferring GWMF from the 3-D temperature distributions. Based on the E2E simulations for gaining zonal mean climatologies of GW momentum flux, we demonstrate that our approach is robust and stable, given a four-track observation geometry and the expected instrument noise under nominal operation conditions. Using phase speed and direction spectra we show also that the properties of individual wave events are recovered when employing four tracks. Finally, we discuss the potential of the proposed observations to address current topics in the GW research. We outline for which investigations ancillary data are required to answer science questions.

1 Introduction

The integration of parameterized gravity waves (GWs) into general circulation models was a tremendous breakthrough in understanding the mesosphere and lower thermosphere (MLT) region. As they replaced Rayleigh friction, the wave-driven circulation and the cold mesopause could be understood (Holton et al., 1995; McIntyre, 1999; McLandress, 1998). However, in their classical formulation, GW parameterizations assume only orography and (often unspecified) non-orographic sources in the troposphere and simplify propagation to be only vertical and instantaneous. In this framework, the waves propagate until they reach either saturation or a critical level and then transfer the dissipated momentum to the background flow only. Such interpretations are supported, for instance, by airglow observations and their match with lower-level filtering as described by, for instance, blocking diagrams (Taylor et al., 1993). In the last 2 decades, however, it has become evident that this view is too simplified and that the simplifications have important consequences for the large-scale dynamics. In order to illustrate this, let us consider three prominent examples in which new concepts are essential: (1) the wind reversal above the summer MLT, (2) the recovery phase of sudden stratospheric warmings, and (3) gravity waves in the thermosphere.

1. At summer midlatitudes tropospheric winds are westerly, but stratospheric winds are easterly. According to the classical picture this should filter out all GWs with phase speeds up to several tens of meters per second in both eastward and westward propagation directions. However, the wind reversal towards westerly winds in the MLT is caused by the dissipation of eastward-propagating GWs. How can these then reach the MLT? One conceivable process would be GWs of extremely high phase speeds (on the order of 90 m s⁻¹) and very small amplitudes, which would not be visible in the stratosphere but would gain saturation amplitudes at high altitude. This occurs in all GW parameterizations with a wide range of phase speeds (e.g., Alexander and Dunkerton, 1999) and has also been suggested by ray-tracing simulations (Preusse et al., 2009a). A second possibility is lateral GW propagation. Indication for lateral propagation of GWs from subtropical convective regions was first found by Jiang et al. (2004) in Microwave Limb Sounder (MLS) observations of the stratosphere. It is seen from the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) observations (see Fig. 5 of Chen et al., 2019) that convective GWs remain at all altitudes in an easterly or low-wind-velocity flow and circumvent the critical levels. Gravity-wave-allowing high-resolution general circulation model (GCM) simulations (H. L. Liu et al., 2014) are consistent with these observations, and the oblique propagation may be favored by large horizontal wind gradients (Thurairajah et al., 2020). The third alternative is secondary (and higher-order) gravity waves. When the original gravity waves from the troposphere break, they exert a body force that excites new gravity waves (Vadas and Fritts, 2002). The relevance of secondary wave generation for the summer MLT is demonstrated by high-resolution GCM simulations (e.g., Becker and Vadas, 2020). There are hence three competing pathways for GWs to reach the summer MLT which need to be distinguished. The way GWs reach the summer MLT necessarily impacts the phase speed and direction distribution: only very fast waves in the first case, a noted poleward preference in the second case, and waves from breaking regions in the third case. It is likely that all pathways occur simultaneously, but the interaction is complex and not well understood (Thurairajah et al., 2020).

2. Gravity waves are believed to play an important role in sudden stratospheric warmings (SSWs) (Thurairajah et al., 2014; Ern et al., 2016; Thurairajah and Cullens, 2022). An SSW event is marked by a general breakdown of the polar vortex and a major event defined by the wind reversal at 10 hPa. A new vortex and an elevated stratopause then form at MLT heights and propagate downward. This re-formation at high altitudes is believed to be mainly caused by GWs. Again the question is how the GWs reach the MLT (Thurairajah et al., 2014; Ern et al., 2016). Here lateral GW propagation and excitation of GWs in the stratosphere, e.g., by an unstable polar vortex, also play a role. Secondary GWs will form at critical levels. For SSWs these additional wave sources not captured by classical GW parameterization schemes jumble the large-scale dynamics in the MLT. In models with parameterized GWs a strong easterly wind bias forms in the MLT after an SSW, which takes the form of a spurious anticyclone (Harvey et al., 2022b) and lasts for more than a month. These results underpin the fact that secondary GWs and other middle-atmosphere sources are essential for the residual circulation and the related dynamical structure in the winter upper mesosphere (Becker and Vadas, 2018; Stober et al., 2021).

3. GWs vertically couple the thermosphere to the lower atmosphere, and thus understanding the wave sources at the lower boundary of the thermosphere is essential for the whole-atmosphere system (e.g., Miyoshi et al., 2014; Park et al., 2014; Yigit and Medvedev, 2015; Yiǧit et al., 2016; Vadas and Becker, 2019; Becker and Vadas, 2020). Wind reversals below the considered altitude are now almost ubiquitous. Still, observations indicate spatial patterns in the global distributions which correlate with those in the stratosphere (Trinh et al., 2018). This is evidence that at least larger parts of the GWs in the thermosphere are secondary GWs, which preserve the spatial patterns of the primary waves. Primary and higher-order GWs reaching altitudes above about 250 km can lead to disturbances or irregularities in the ionospheric layer and thereby affect space-based applications (e.g., Hines, 1960; Bertin et al., 1975; Vadas and Fritts, 2006; Vadas, 2007; Krall et al., 2013; Nishioka et al., 2013; Yiǧit et al., 2016; Liu, 2016). Together with tides and planetary waves, GWs are the most important dynamical process in the MLT region. Understanding the various aspects regarding GW instability and transition turbulence, interactions with the ambient large-scale flow, and the generation of higher-order GWs requires extensive knowledge of the spectral and spatial distributions of GWs, including geographical and seasonal variations.

This leads us to some higher-level science questions; answering them is essential for understanding the MLT and the coupling between the middle atmosphere and the thermosphere.

Science questions are the following.

• How do GWs reach the summer MLT?

• Which GWs lead to the formation of the elevated stratopause and the new vortex after an SSW?

• Is there strong westward GW drag in the MLT in the period 2 weeks to 2 months after an SSW?

• Which GWs propagate into the MLT?

In order to answer these questions we need to characterize the GWs in the MLT region. Key quantities to be determined are

• zonal mean GW momentum flux and its vertical gradient (GW drag) as well as

• phase speed and direction distributions of GWMF.

In order to close the momentum budget, the zonal mean of the zonal GW momentum flux is particularly required, but zonal mean meridional momentum flux may contribute as well (Ern et al., 2013a). In order to calculate a meaningful zonal mean, global coverage is required. Phase speed spectra are an essential tool to quantify the interaction with the background wind (e.g., Taylor et al., 1993). The spectra are representative for an evaluation region and can be constructed from individual GW observations. Still, for the global picture global coverage is required as well. This can be only provided from satellite observations.

For our study the zonal mean of zonal GW momentum flux is of particular importance as the values directly inferred from the winds provide a true reference value. This is, to a somewhat lesser degree, also true for the meridional momentum flux, as will be discussed below.

In general, no observation technique can characterize the entire spectrum of GWs, and different kinds of observations need to be combined for a consistent picture of GWs and their impact in the MLT. For a limited number of locations, spectral information and GW momentum flux can be inferred from ground-based radar and lidar systems (e.g., Stober et al., 2013; de Wit et al., 2014; Placke et al., 2015; Bossert et al., 2015, 2018; Chum et al., 2021). In addition, ground-based airglow imagers provide information about GWs with long vertical and short horizontal wavelengths (e.g., Tang et al., 2002; Espy et al., 2006; Shiokawa et al., 2009). When it comes to the large-scale momentum budget in the MLT, these observations are, however, biased. They are made only on land and often in locations of specific geophysical interest (for example, strong activity of mountain waves). Furthermore, optical systems can work under clear-sky conditions only.

Figure 1 Proposed observation geometry of two CubeSats flying in parallel and viewing backward. Each CubeSat carries an SHI which images the atmosphere in the vertical and generates a combined spectral–spatial view in the horizontal. By splitting the interferogram at zero optical path difference, both sides can be evaluated individually, allowing for four effective observation tracks.

Based on existing spaceborne limb-scanning observations that allowed distributions of the absolute GW momentum flux to be inferred (e.g., Ern et al., 2004; Alexander et al., 2008; Preusse et al., 2009a; Ern et al., 2011; Alexander, 2015; Ern et al., 2018), proposals were made on how a limb-imaging satellite mission could drastically improve our knowledge about GWs in the stratosphere (Preusse et al., 2009b, 2014). Such an instrument would provide 3-D data at good spatial resolution by high along-track sampling, tomographic retrieval in along-track slices (Ungermann et al., 2010a, b; Song et al., 2017), and across-track coverage by multiple tracks (illustrated in Fig. 1). Still, general restrictions due to the radiative transfer along the line of sight remain despite tomographic retrievals, and the observational filter allows observing only waves with horizontal wavelengths longer than 100 km (Preusse et al., 2009b). The first existing global observations exploited for 3-D data were nadir measurements of the Atmospheric InfraRed Sounder (AIRS) (Ern et al., 2017; Hindley et al., 2020). Although these data capture only long vertical wavelengths, i.e., high intrinsic phase speeds, they provide information on direction characteristics and allow demonstrating how backward ray tracing can be used for source identification from global data (Perrett et al., 2021). These examples are for the stratosphere only. Nevertheless, it is evident from such studies that a limb imager would provide novel information about GWs in the MLT.

The spatial sampling drives the complexity of the instrument and hence drives the cost of the mission. The spatial sampling requirements therefore need to be justified. The across-track dimension is provided by observing multiple tracks. An important question is therefore how many parallel tracks are required to gain reliable information about medium-scale GWs with horizontal wavelengths longer than 100 km (i.e., the ones visible to a limb sounder) and how these tracks should be spaced. On the one hand, it is obvious that the wider the overall swath is and the smaller the individual pixels are, the higher the likelihood is to acquire unprecedented scientific data. On the other hand, a larger number of tracks is a driver for increased instrument complexity and data downlink capacity. Besides the traditional satellite missions with a cost easily on the order of several tens of millions of Euro, an alternative option would be a CubeSat mission which takes fewer tracks of measurements but is still capable of providing a similar amount of information about GWs. Employing a spatial heterodyne spectrometer (SHS), the CubeSat instrument will provide good spectral resolution in the selected emission band, which is helpful to constrain the retrieval and obtain accurate temperature. Compared to an imager, however, the number of parallel tracks is lower. One of the aims of this paper is therefore to examine the minimum number of tracks required for deriving 3-D wave vectors of GWs from tomographic temperature observations.

From the above motivation, we deduce the following observation concept: airglow emissions at 762 nm from the oxygen atmospheric band (O₂ A band) are particularly suited to gain information on MLT dynamics. Limb observations facilitate high vertical resolution during both day and night. Assuming rotational local thermodynamic equilibrium, the kinetic temperature around the tangent points can then be inferred from the relative line intensities. Using advances in CubeSat standard components, detector technology, and optics, a highly miniaturized spatial heterodyne interferometer (SHI) (Kaufmann et al., 2018) was developed for this purpose. This detection technology can be applied in a CubeSat constellation mission, which consists of two CubeSats, each hosting an SHI. By flying the two SHIs in parallel (illustrated in Fig. 1) and by splitting one interferogram into two left-hand and right-hand parts and separately mirroring each parts (see Sect. 3.5), thus splitting the horizontal field of view (FOV) in two, in total four independent observation tracks can be obtained from the proposed satellite observation geometry. High along-track resolution will be achieved using tomographic retrievals.

Based on limb-sounding four tracks simultaneously, we aim to observe medium-scale GWs of horizontal wavelengths longer than 100 km. In order to demonstrate the feasibility of quantifying GW properties by airglow limb observations, we will perform a full chain of end-to-end (E2E) simulations based on self-consistent model data. This validation of our methodology is based on GW parameters (e.g., GW momentum flux) that are derived from temperature residuals using the GW polarization relations (Fritts and Alexander, 2003; Ern et al., 2004; Preusse et al., 2009b). The feasibility will then be evaluated by inferring such GW parameters directly from the model data. This means that we choose the inferred GW parameters as a performance measure instead of a separate consideration of noise and resolution¹. The stringency of validation by comparing distributions of GW momentum flux with reference distributions is discussed in Sect. 2.3. Phase speed spectra are essential for understanding the interaction with the background winds. However, all spectral investigations are dependent on the choice of the method and hence no absolute reference exists. We therefore compare in this paper how the spectra degrade when fewer tracks are employed.

In order to study the viability of a CubeSat mission, we will consider the following questions in this study.

• Are polarization relations valid in the MLT region?

• How few measurement tracks are required? Are four or even two tracks sufficient?

• How much instrument noise can we afford in the temperature retrieval and wave analysis?

We address these three questions as follows: the assessment strategy is outlined in Sect. 2. Detailed introductions to models, tools, and the instrument are presented in Sect. 3. The outcomes of the assessment and answers to the questions are given in Sect. 4, followed by a discussion on scientific applications in Sect. 5. Finally, we summarize our findings (Sect. 6).

2 Assessment strategy

In this section we outline the strategy to address the questions formulated in the Introduction. We base our study on fields of winds and temperature from a free-running general circulation model (GCM) called the HIgh Altitude Mechanistic general Circulation Model (HIAMCM) (Sect. 3.1) resolving a larger part of the GW spectrum. Since the GCM simulates all dynamical features from first principles, wind and temperature structures of the model fields are consistent for all waves resolved. This consistency is essential for diagnosing deviations between the reference, analyzed GW momentum flux, and full end-to-end simulations: the smaller the deviations, the better the performance of the observation system. It should also be noted that while a realistic representation of the actual atmosphere by the model is important, potential deviations of the model fields from the simulated meteorological situation are not part of the assessment. To retain this consistency, winds and temperatures are separated into global-scale dynamics and GW fluctuations in the same way on the full model fields (Sect. 2.1). For an assessment we need a reference to evaluate the performance against. This is the zonal mean of zonal momentum flux derived from the model winds (Sect. 2.2). With this given, we can design the method (Sect. 2.3) to tackle the three questions posed in the Introduction.

2.1 Scale separation

The assessment is based on consistent GW-related fluctuations of temperature and wind velocities. This requires a scale separation between GWs on the one hand and global-scale dynamics on the other hand. This is usually performed by separating by zonal wavenumber.

For satellite data, traditionally waves up to zonal wavenumber 6 have been treated as global-scale waves and all remaining fluctuations as GWs (e.g., Fetzer and Gille, 1994; Preusse et al., 2002; Ern et al., 2018). This separation is used as wavenumber 6 is the highest wavenumber which can be reliably resolved by a single-observation-track low Earth orbit (LEO) satellite (Salby, 1982). Such a detrending method via space–time spectral analysis has been used for, e.g., SABER data (Ern et al., 2018).

For model studies a much higher separation wavenumber has often been used. Strube et al. (2020) have shown that separation wavenumbers 6 to 8 are sufficient for the stratosphere and that removing wavenumbers up to 40 significantly cuts into the GW part of the spectrum. For studies including the UTLS Strube et al. (2020) recommend zonal wavenumber up to 18 for scale separation. In this study we also use zonal wavenumber 18 with an additional meridional Savitzky–Golay filter of third-order polynomials over 5^∘ of latitude. This defines our large-scale background, which is subtracted from individual temperature values in order to define residuals.

A satellite in a sun-synchronous orbit acquires data at a continuously evolving observation time but fixed local time for a given latitude on ascending and descending orbit legs. Such a sampling cannot be generated from model data which are sampled at fixed UTC and sampling intervals of O(1 h). Switching between model fields as the orbit evolves would result in jumps at the switching points, while interpolation in time would smooth the interpolated GW fields in an unpredictable manner. In this study we therefore use a single model snapshot (1 January 2016 at 06:00 UT). Synthetic orbit data generated in this manner allow addressing all questions stated in the Introduction. In particular, a fixed UTC has the advantage that the synthetic orbit data can also be compared to the reference of zonal mean GWMF from full model fields, which is the basis of our assessment.

2.2 Zonal mean momentum flux: a true reference

At the end, we aim to quantify the vertical flux of horizontal pseudomomentum of GWs (Fritts and Alexander, 2003):

$$\begin{array}{}\text{(1)}& \left({\mathit{F}}_{px},{\mathit{F}}_{py}\right)=\mathit{\rho }\left(\mathrm{1}-{\displaystyle \frac{f^{\mathrm{2}}}{{\widehat{\mathit{\omega }}}^{\mathrm{2}}}}\right)\left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{\mathit{w}}^{\prime }},\stackrel{\mathrm{‾}}{{\mathit{v}}^{\prime }{\mathit{w}}^{\prime }}\right),\end{array}$$

where F_px indicates its zonal component and F_py is its meridional component. ρ is the background density, f the Coriolis parameter, $\widehat{\mathit{\omega }}$ the intrinsic frequency, and u^′, v^′, and w^′ the wind vector perturbations due to the GW in zonal, meridional, and vertical directions, respectively; eastward, northward, and upward are positive signed. The overline denotes the average over a full or multiple wavelengths of the wave. The vertical gradient of the pseudomomentum flux (PGWMF) determines the acceleration of the background wind on a rotating sphere. However, determining $\widehat{\mathit{\omega }}$ from a given 3-D data set involves some kind of wave analysis and, accordingly, assumptions. On the other hand, the zonal mean of zonal momentum flux ($\mathit{\rho }\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{\mathit{w}}^{\prime }}$) is a true reference as it depends on the wind fluctuations only and as on the cyclical domain of longitude all waves are properly averaged.² If not explicitly stated otherwise, we will hence consider GW momentum flux (GWMF) without the correction for Coriolis force in our assessment.

2.3 Method of assessment

The E2E approach for evaluating the performance of the proposed mission concept is illustrated in Fig. 2. We employ HIAMCM data (Becker and Vadas, 2020, see Sect. 3.1) for a realistic, self-consistent basis of the E2E simulations. The distributions of temperature, density, and wind velocities are separated for large-scale structures such as planetary waves and tides on the one hand and small-scale residuals due to GWs on the other hand. We consider both wind and temperature data. Temperatures are the observation target of the measurement method. From the winds we gain our reference for the assessment. As described in the Introduction, zonal mean GWMF calculated directly from the wind perturbations provides an unambiguous reference of truth against which we can compare the values from simulated observations and thus quantify the influence of the various assumptions or constraints needed for observing GWMF with a real instrument.

Figure 2 Schematic flow diagram of the assessment steps to address the three major questions highlighted in yellow diamond boxes. The assessment is based on comparison with zonal mean GWMF calculated directly from the model winds and considered here as the reference of truth. From the top down further assumptions and/or constraints are added as the tested data become more similar to the real observations.

The first question which we address is the applicability of the polarization relation for the MLT. GWMF can be deduced from 3-D temperature data by determining the 3-D wave vector assuming polarization and dispersion relations (see Sect. 4.1). In order to test the validity of this approach, we apply our wave analysis – the small-volume 3-D sinusoidal fit method (S3D; Sect. 3.7), which provides wave amplitudes and 3-D wave vectors based on small sub-volumes of the data set considered – directly to the full model data and compare thus generated zonal mean GWMF with our reference.

The second question regards the number of tracks required for reliable GWMF quantification (see Sect. 4.2). For this assessment, we sample model residual temperatures onto various synthetic measurement geometries for different numbers of observation tracks, apply S3D and, again compare against the reference.

Third, we assess the influence of the observation technique which comprises both the observational filter as a general limitation of limb sounding and the specific instrument parameters (see Sect. 4.3). This assessment therefore encompasses all steps of observation starting from synthetic radiances generated by forward modeling, a simplified instrument model adding realistic noise, and a tomographic retrieval. Again, S3D is applied to the outcome of the full E2E simulations, and zonal mean GWMF is assessed against the reference. Different noise levels are tested in order to determine the robustness against instrument performance.

3 Models, tools, and instrument

This section describes the atmospheric circulation model and the radiative transfer model we base our study on, introduces the SHS instrument proposed for the observation, and gives an introduction to the wave analysis and retrieval tools employed.

3.1 HIAMCM

The HIAMCM (HIgh Altitude Mechanistic general Circulation Model) (Becker and Vadas, 2020; Becker et al., 2022) is a new high-altitude version of the KMCM (Kühlungsborn Mechanistic general Circulation Model). The HIAMCM employs a spectral dynamical core with a terrain-following hybrid vertical coordinate. This includes a correction for non-hydrostatic dynamics and physically consistent thermodynamics in the thermosphere. The current model version uses a triangular spectral truncation at a total wavenumber of 256, which corresponds to a horizontal grid spacing of ∼ 52 km. The altitude-dependent vertical resolution includes 280 full levels. The level spacing is dz ∼ 600 m below z ∼ 130 km, with dz increasing with altitude farther above and dz ∼ 5 km above z ∼ 300 km.

The HIAMCM captures atmospheric dynamics from the surface to approximately 450 km. GWs are simulated explicitly with an effective resolution that corresponds to a horizontal wavelength of ∼ 200 km. Non-resolved scales are parameterized by macro-turbulent vertical and horizontal diffusion based on the Smagorinsky model. Since molecular viscosity is taken into account for both vertical and horizontal diffusion, the HIAMCM does not require an artificial sponge layer. Resolved GWs are dissipated self-consistently by molecular diffusion in the thermosphere above ∼ 200 km and predominantly by macro-turbulent diffusion at lower altitudes. These features allow the HIAMCM to capture the generation, propagation, and dissipation of medium-scale GWs, including their interactions with the large-scale flow and the generation of secondary and tertiary waves. This capability is essential to simulate GW dynamics in the MLT (Becker and Vadas, 2018) and at higher altitudes (Vadas and Becker, 2019).

The HIAMCM employs radiation and moist convection schemes that are simplified compared to comprehensive methods. Though convection is only present in the troposphere, it is essential as a source of upward-propagating waves, namely GWs, tides, and tropical wave modes, and is thus also highly relevant for the dynamics in the MLT. Furthermore, the model does not include a chemistry module, and ion drag is the only ionospheric process that is accounted for. To distinguish these idealizations from methods employed in community models, the HIAMCM is said to be a mechanistic model. Nevertheless, the key features of a climate model (topography, simple ocean model, radiative transfer, boundary layer processes, tropospheric moisture cycle) are fully taken into account. Also note that the HIAMCM is currently the only GW-resolving whole-atmosphere model that can be nudged to reanalysis in the troposphere and stratosphere, allowing for the simulation of observed events (Becker et al., 2022).

3.2 Orbit simulator

The simulation of the orbit (illustrated in Fig. 3) is based on a fixed orbit inclination, orbit altitude (shown in Table 1), and start longitude at the beginning of the day. Assuming a spherical Earth and constant gravity acceleration scaled to orbit altitude the in-orbit velocity is determined. A time series of orbit positions on Earth surface is calculated from the satellite position on the orbit fixed in space and the rotation of the Earth. A grid for atmosphere representation (atm-grid in the following) is then generated based on the tangent points of the backward-viewing direction for 80 km altitude and spanning a local rectilinear grid with the x direction along this tangent point track, a y direction perpendicular to this tangent point track, and the local vertical. Both satellite position and atm-grid are then used to build a complete set of matching tangent points for radiative transfer and retrieval simulations as described in Sect. 3.3 and 3.6. The basis for the radiative transfer simulations is the atmospheric quantities such as temperature and pressure interpolated from the longitude–latitude grid of the general circulation model to the atm-grid by means of spline interpolation.

Figure 3 Schematic view of the synthetic orbit and grid generation.

Table 1 Orbit parameters used for the simulated observations

3.3 O₂ A-band airglow photochemistry and radiative transfer process

This section describes the photochemical processes which initiate the generation of the O₂ A-band emission, followed by a short description of radiative transfer process propagating the radiance to the instrument.

Figure 4 (a) Schematic of the production of excited-state ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$; dashed lines indicate neglected transitions. (b) Number density of excited O₂ molecules due to the different production mechanisms using HAMMONIA model data (Schmidt et al., 2006); note that during daytime all five production mechanisms are active, whereas during nighttime only the Barth process is active.

O₂ can be present as ¹⁶O₂, ¹⁷O¹⁶O, or ¹⁸O¹⁶O; the latter two can be neglected due to their low abundance, as shown by Slanger et al. (1997). In general, an excited molecule can be in one of multiple electronic states; for each electronic state the radicals can be in one of multiple vibrational states. The transitions between different electronic–vibrational states form atmospheric emission and absorption bands. Each band consists of multiple emission lines due to the transitions within multiple rotational states of vibrational states. We measure the A-band emission at 762 nm, which is the electronic transition from the second excited state ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$ to the ground state ${\mathrm{O}}_{\mathrm{2}}(X^{\mathrm{3}}{\mathrm{\Sigma }}_g^{-},v=\mathrm{0})$. A detailed description of the dayglow O₂ A band is given by Sheese (2009), Bucholtz et al. (1986), and Zarboo et al. (2018). It has three sources, which are depicted in Fig. 4a.

First, the excited state can be produced by photon absorption in the atmospheric bands. Bucholtz et al. (1986) show that the γ-band absorption from ${\mathrm{O}}_{\mathrm{2}}(X^{\mathrm{3}}{\mathrm{\Sigma }}_g^{-},v=\mathrm{0})$ to ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{2})$ can be neglected, and hence it is not shown in Fig. 4a. Only the absorption in the A and B bands is therefore considered. The excited molecules in ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{1})$ are rapidly deactivated to ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$ via a quenching process. Slanger et al. (1997) argue that the B-band emission at 686 nm is insignificant compared to the quenching process and can be neglected. The emission in the (1,1) band is considered. The second source is due to photolysis of O₂ in the Schumann–Runge continuum (J_SRC) and at Lyman α ($J_{{\mathrm{L}}_{\mathit{\alpha }}}$) and due to the photolysis of O₃ in the Hartley band (J_H). It produces excited atomic oxygen O(¹D), which transforms to ${\mathrm{O}}_{\mathrm{2}}\left(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+}\right)$ due to collisional excitation with O₂ in ground state. The third source is a chemical source, which is independent of solar radiation and hence also present during nighttime. The source of this process is a three-body recombination of atomic oxygen, producing electronically exited ${\mathrm{O}}_{\mathrm{2}}^{\ast }$ radicals. From there, ${\mathrm{O}}_{\mathrm{2}}\left(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+}\right)$ is produced through a direct quenching process or a chain of quenching processes, going through O(¹S) and O(¹D). This process was first described by Barth and Hildebrandt (1961) and is called the Barth process. Since some of the related rate coefficients are not well known, McDade et al. (1986) proposed a model with fitting parameter to describe the Barth process. A detailed description of the calculation of O₂ A-band emission is given in Appendix A. Figure 4b shows the number density of excited O₂ molecules due to the different production mechanisms using HAMMONIA model data (Schmidt et al., 2006). Note that during daytime all five production mechanisms are active, whereas during nighttime only the Barth process is active.

The O₂ A-band airglow emissions are transmitted through the atmosphere before they are detected by an instrument. The observed airglow spectra are integrated slant-path radiances along the instrument viewing line of sight (LOS). For limb observations, the emissions from the lowermost tangent points along the LOS contribute the most to the integrated radiances provided that the atmosphere is still optically thin for the emission lines.

Due to the high abundance of the ground-state O₂ molecule in the atmosphere and its self-absorption effect, the emitted O₂ A-band radiance can only partly pass through the atmosphere and reach the instrument detector. This radiative transfer process is described by the Lambert–Beer law, and the observed spectral irradiance intensity can be written in integral form as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}I\left(\mathit{\nu }\right)=& \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}I\left(s\right)D(\mathit{\nu },s)\\ & \exp \left(-\underset{s}{\overset{\mathrm{\infty }}{\int }}n\left(s^{\prime }\right)\mathit{\sigma }\left(s^{\prime }\right)D(\mathit{\nu },s^{\prime })\mathrm{d}{s}^{\prime }\right)\mathrm{d}s,\end{array}\end{array}$$

where ν refers to the wavenumber of the spectral line, s denotes the LOS distance, D(ν) represents the line shape broadening profile, which is dominated by Doppler broadening in the middle and upper atmosphere, n(s) is the number density of O₂, and σ is the O₂ absorption cross-section.

In the lower atmosphere, nearly all O₂ A-band emissions are self-absorbed. Thus, O₂ A-band airglow cannot be detected on the ground. During daytime, the O₂ A-band emissions are sufficiently strong for observations at altitudes from 60 to 120 km and are reduced to the range of 80 to 100 km during nighttime. In the simulation, we have assumed optically thick conditions for the center wavelength of the emission lines (see Appendix B).

3.4 Spatial heterodyne interferometer

A spatial heterodyne interferometer is based on the principle of a Michelson interferometer, but the two mirrors are replaced by fixed tilted gratings. The simplified concept of the SHS is shown in Fig. 5. Light of a frequency within the spectral bandpass enters the instrument along the optical axis. Light of each frequency is split into two waves by the beam splitter and diffracted and reflected by the gratings. The two waves arrive back to the beam splitter, producing an interference pattern, which is forwarded to the detector. Considering multiple emissions within the bandpass, an interferogram consists of multiple superposed interference patterns along the x axis, which can be described by sinusoidal waves. The spatial frequency of the produced fringe pattern follows from the grating equation, denoted by

$$\begin{array}{}\text{(3)}& \mathit{\sigma }(\sin {\mathit{\theta }}_{\mathrm{L}}+\sin ({\mathit{\theta }}_{\mathrm{L}}-\mathit{\gamma }\left)\right)={\displaystyle \frac{m}{d}},\end{array}$$

where σ is the wavenumber of the incoming light, θ_L the blaze angle of the gratings, also called the Littrow angle, d⁻¹ is the grating groove density, m is the diffraction order, and γ the outgoing diffraction angle with respect to the Littrow angle. The Littrow wavenumber is the wavenumber where γ=0, thus

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathrm{L}}={\displaystyle \frac{m}{\mathrm{2}\sin \left({\mathit{\theta }}_{\mathrm{L}}\right)d}}.\end{array}$$

Figure 5 Schematic of the SHS instrument.

Using Taylor expansion, Harlander et al. (1992) show that the spatial frequency is dependent on the wavenumber by

$$\begin{array}{}\text{(5)}& f\left(\mathit{\sigma }\right)=\mathrm{4}(\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{L}})\tan \left({\mathit{\theta }}_{\mathrm{L}}\right)M,\end{array}$$

where M is the magnification factor introduced by the camera optics.

This equation shows that the relation between spatial frequency and wavenumber is symmetric around the Littrow wavenumber.

Following Roesler and Harlander (1990) and Deiml (2017), an ideal one-dimensional interferogram along the x axis can be described by

$$\begin{array}{}\text{(6)}& I\left(x\right)=\underset{b_{\mathrm{0}}}{\overset{b_{\mathrm{1}}}{\int }}B\left(\mathit{\sigma }\right)\left[\mathrm{1}+\cos \left(\mathrm{2}\mathit{\pi }f\right(\mathit{\sigma }\left)x\right)\right]\mathrm{d}\mathit{\sigma },\end{array}$$

where B is the spectral radiance, and b₀ and b₁ are the lower and upper bound of the spectral filter, respectively.

The interferogram is transformed into a spectrum by Fourier transformation. Figure 6 shows the O₂ A-band emission for two temperatures and the corresponding spectra as seen by the instrument. The temperature dependency as seen in the lower panel of Fig. 6 is used to retrieve temperature.

Figure 6 Modeled O₂ A band seen by the SHI instrument for temperature at 200 and 210 K. The spectra are presented by the Dirac impulses such that the sum over all emission lines is equal to 1; below the relative intensity difference between the spectrum at 200 and 210 K is shown. The dashed line shows the theoretical filter curve and the dotted vertical line the designed Littrow wavenumber of the AtmoLITE instrument. The gray shaded area shows the spectrum convolved by the instrument line shape; r = 1.43 cm⁻¹ is the spectral step size.

A silicon-based detector is used in this instrument. The operation in ambient to cool conditions gives a shot-noise-limited system for integration times of 1–10 s (Liu et al., 2019). Shot noise can be modeled by a Poisson process with mean and variance equal to the signal. For a signal above 10 counts, the Poisson distribution approximates a normal distribution about its mean. Thus, for simplicity the shot noise is approximated by additive white Gaussian noise with standard deviation equal to the square root of the signal in each pixel.

3.5 Split of the interferogram

Revisiting Fig. 5, the front optics map the captured atmospheric scene onto the gratings. In a similar fashion, the camera optics map the image on the grating to the image on the detector. Accordingly, the SHI performs a spatial mapping of the atmospheric scene onto the detector. As described in Sect. 3.4, the spectral information is contained on the x axis due to constructive and destructive interference induced by the gratings. The interferogram therefore contains spatial information on the vertical (y axis) and superimposed spectral and spatial information on the horizontal axis (x axis). This is utilized when splitting the interferogram at zero optical path difference (ZOPD) and using each half individually. The spectral information is fully contained in one half as the interferogram is symmetric around the ZOPD by definition. Mirroring each side around the ZOPD gives a full interferogram for each side, which can be employed to derive two independent temperatures along the horizontal axis. Each interferogram then contains temperature information of the associated side of the field of view. The concept to mirror the interferogram at the ZOPD has already been used by Johnson et al. (1996) for the far-infrared spectrometer (FIRS)-2 and by Gisi et al. (2012) for the TCCON FTIR spectrometer to gain a higher resolution. However, Ben-David and Ifarraguerri (2002) and Brault (1987) point out that the phase correction including finding the correct ZOPD is crucial.

Figure 7 Simulation of interferogram split assuming a linear temperature gradient. Panel (a)​​​​​​​ shows the temperature gradient from 200 to 220 K across the horizontal field (solid line) with each side symmetrically extended around the center (dashed lines); the red circle and the blue diamond indicate the location of the retrieved temperature for the left and right hand side. Panel (b) shows the left symmetrically extended interferogram. Panel (c) shows the temperature retrieval of (b). Panel (d) shows the difference between the left and right symmetrically extended interferogram. Panel (e) shows the difference between the left and right symmetrically extended raw spectrum.

The following simulation demonstrates that averaged temperatures of the respective parts of the field of view can be recovered from an analytical interferogram, theoretically allowing two independent cross-track measurements with one instrument. We assume a linear temperature gradient across the horizontal field. We simulate each pixel with the associated temperature and assemble the full interferogram pixel by pixel. Subsequently, the interferogram is split at the ZOPD, and each side is symmetrically extended to get two full interferograms. Each side can be then used to retrieve a temperature. Figure 7 shows the simulation result of an ideal interferogram without noise using a simple temperature retrieval which minimizes the squares of the residuals. The red circle and the blue diamond in Fig. 7a indicate the retrieved temperature. Mapping the retrieved temperatures back onto the initial temperature gradient results in two data points which are about 17 km apart. Note that the two retrieved temperatures are shifted towards the center and are not located at the spatial centers of each side at ±15 km. The Fourier transformation is a weighted sum over the samples in the interferogram. Areas which deviate the most from the mean of the interferogram therefore contribute the most to the shape of the spectrum. Since the area around the ZOPD deviates the most from the interferogram's mean, the temperature information of that area contributes the most to the retrieved temperature. This causes the shift towards the center. An in-depth validation using calibration and orbit data considering noise and instrument errors is the content of future research.

3.6 Tomographic retrieval for generating 3-D atmospheric volumes

Retrieving temperatures from measured limb spectra is a classic inverse problem. That is, we have a radiative transfer model that can compute measured spectra from an assumed atmospheric state (forward model), but the inverse problem is much harder, since it is typically both underdetermined (multiple atmospheric states could result in the same set of measurements) and ill-posed (there might be, in theory, no atmospheric states that could result in a given real-life measurement affected by instrument noise and other error sources). Such a problem is solved using the forward model and a mathematical framework for inverse modeling (e.g., Rodgers, 2000). The main idea of this approach is iterative minimization of the following function (called cost function):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}J\left(\mathit{x}\right)& ={\left(\mathit{y}-\mathbf{F}\left(\mathit{x}\right)\right)}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\left(\mathit{y}-\mathbf{F}\left(\mathit{x}\right)\right)\\ & +{\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right)}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{apr}}\right).\end{array}\end{array}$$

Here y is a set of measurements taken by the instrument, and x is the candidate atmospheric state. F(x) is the forward model; it maps an atmospheric state x to the set of measurements that the instrument would acquire if the atmosphere was indeed in the state x. x_apr represents our prior knowledge of the atmospheric state. In our case that is an estimate of the large-scale temperature structure of the atmosphere without gravity waves. ${\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}$ and ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ are positive-definite symmetric matrices (covariance matrices).

The first term on the right-hand side of Eq (7) quantifies how closely the atmospheric state x matches the observations y. One can construct the covariance matrix ${\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}$ knowing the measurement error characteristics typical for the instrument in question. The second term quantifies how likely the atmospheric state x is given our prior knowledge about the atmosphere. This knowledge includes both the base state of the atmosphere x_apr and more general considerations (such as the fact that large spatial discontinuities in temperature are unlikely) that are taken into account when constructing ${\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}$ (Ungermann et al., 2010a). Equation (7) is solved iteratively using the Levenberg–Marquardt algorithm (Marquardt, 1963) and a conjugate gradients solver. For both the forward model and the inverse modeling we need a computationally efficient implementation.

Here, we employ the JURASSIC2 forward model (e.g., Ungermann, 2013) to simulate the spectra based on a 2-D discretization of the atmosphere along the satellite track. This model was enhanced for this study by a simple adjoint line-by-line model dedicated to the simulation of the O₂ A band (see Sect. 3.3). The inversion uses the JUTIL Python library to bring the synthetic and measured spectra in agreement using a truncated conjugate gradient trust region method (e.g., Ungermann et al., 2015).

In order to perform a full end-to-end test, noise is added to simulated observations to emulate instrument performance. To generate the noise in the synthetic measured spectra, we compute from the spectrum the number of photons hitting (on average) one detector pixel per second and use this number to determine Gaussian noise assuming that the instrument is shot-noise-limited and the dark current is negligible; i.e., the 1σ is computed from the square root of number of photons. Note that we assume an efficiency of 0.2; i.e., only one-fifth of inbound photons will end up in the modulated part of the interferogram. This noise is then reduced according to assumed averaging in time and space (4.5 s integration time, 13 vertical detector rows, four horizontal detector columns).

The individual measurement tracks of the proposed satellites are treated separately from one another, and each constitutes a separate 2-D slice at the orbital plane of the satellite or very close to it (Fig. 8). Limb measurement geometry and backward-viewing direction result in each air parcel within the 2-D slice being observed from multiple points of the satellite orbit. This allows reconstructing the 2-D temperature cross-section in a tomographic fashion. The full 3-D state is reassembled afterwards from the individually retrieved 2-D slices. The satellite speed allows gathering all relevant measurements for a spatial sample in the order of minutes, which is short compared to typical periods of gravity waves observable by our instrument.

Figure 8 The principle of 2-D tomography. The field of view of the satellite within the orbital plane is shown for different satellite positions along the orbit. A point A inside the airglow layer lies within the field of view for all three positions, and therefore airglow radiance at point A contributes to observed radiances at all these positions and every position in between. We can use the observations to solve for airglow radiance at point A (and all other such points) using inverse modeling.

As explained in Sect. 3.3 there are more excitation mechanisms and higher emissions during daytime than during nighttime. For one retrieval, we therefore split the orbit at the position of the terminator in roughly two equal halves, one containing the daytime and the other containing the nighttime measurements (even though we currently apply the same retrieval settings to both). We then lengthen the two orbit parts by an additional ≈ 1500 km on both ends. Horizontally, a sampling of 30 km is used, while the vertical range is sampled in 1 km steps from 60 to 120 km and in 2 km steps above (Table 2).

Table 2 Summary of detector and sampling properties.

^∗ The spectral wavenumber is defined here as λ⁻¹, where λ is the wavelength of emissions or absorptions.

The forward simulations employ a regular grid of two rays per detector row, which are combined assuming a line spread function of Gaussian shape with 1σ corresponding to the height of the row.

The spectra range from 13 060 to 13 160 cm⁻¹ with a sampling distance of 2 cm⁻¹ and a spectral resolution (full width at half-maximum, FWHM) of ≈ 3.9 cm⁻¹ due to the employed strong Norton–Beer apodization (Norton and Beer, 1976). More details about the simplifying assumptions for the forward model in the tomographic retrieval are given in Appendix B.

3.7 JUWAVE S3D wave analysis

In this study we investigate GWs in a narrow stripe (small swath width) of observations along the tangent point tracks. This requires an analysis method which can, at least in one direction, analyze waves with notably longer wavelengths than the size of the analysis volume. In addition, the vertical wavelength of gravity waves is refracted by the background wind; this contradicts the assumption of a stationary wave spectrum over the range of the analysis volume, which is made, e.g., by Fourier transform. For such applications, the small-volume sinusoidal fit method (S3D) was developed and tested for the purpose of GW analysis in small observation volumes and for highly localized GW fields (Lehmann et al., 2012).

In this method, the observation volume or model domain is divided into small sub-volumes and a sinusoidal fit is performed on each sub-volume:

$$\begin{array}{}\text{(8)}& T_i^{\prime }=\sum _j{A}_j\sin \left({\mathit{k}}_j{\mathit{x}}_i\right)+B_j\cos \left({\mathit{k}}_j{\mathit{x}}_i\right),\end{array}$$

where $T_i^{\prime }$ is the temperature fluctuation at the location x_i, and k_jx_i is the scalar product between the wave vector of the j wave component with the spatial coordinate vector of the i point in the analysis volume. The wave components j are determined sequentially by subtracting the wave field of component j from the fit volume before fitting j+1. In this study, three wave components are fitted. Amplitudes and wave vectors are determined via least squares fit: amplitudes are determined analytically and wave vectors via a variational approach. The minimum χ² from a steepest descent method and a nested interval method is selected.

The size of the analysis cube is selected in a way that most of the spectral content has a wavelength of ${\mathrm{cs}}_{\mathrm{d}}/\mathrm{2}<{\mathit{\lambda }}_{\mathrm{d}}<\mathrm{3}{\mathrm{cs}}_{\mathrm{d}}$, with λ_d being either the horizontal or the vertical wavelength and cs_d the analysis volume diameter. This choice is motivated by previous sensitivity studies (Preusse et al., 2012) and will be discussed further in Sect. 4.2.1. In particular, we find that a vertical cube size of 15 km comprises the spectral power in the MLT region almost entirely. In the lower thermosphere, which is analyzed for consistency reasons as well, larger cube sizes are needed. To retain only reliable fits, we omit all fits with derived horizontal or vertical wavelengths larger than 3 times the respective cube size from evaluation. In order to enhance the vertical resolution and hence to better capture the loss of GWMF by the approach to a critical level at the MLT wind reversal, a refit of only amplitude and phase based on the wave vector from the initial fit is performed. For this study, we keep the horizontal cube size the same but reduce the vertical cube size to 5 km, with the exception of Sect. 4.1, where the results are obtained through a vertical cube size reduction to 4 km. In the same section, the initial cube grid consists of cubes with sizes 300 km × 300 km × 15 km (at 75 km altitude) and 600 km × 600 km × 20 km (at 130 km altitude).

The synthetic observation data have a fixed sampling in the x, y, and z direction, on which the analysis cube size is defined via the number of sampling points. For the model data, a fixed model sampling in terms of degrees longitude in the zonal direction means a coarser (in distance) sampling close to the Equator and a finer sampling at high latitudes due to the shorter distance between two respective longitudes at higher latitudes. Therefore, the size of a fixed cube is specified in kilometers instead of degrees and the number of fitting points is adapted accordingly. This ensures that the same part of the spectrum is targeted independent of latitude along the longitude direction.

The results of S3D are expected to be a good compromise between spatial and spectral representation. It has been shown by Lehmann et al. (2012) that the spectral distribution composed of all S3D wave fits in a given region reproduces the spectral content obtained via Fourier analysis of the same region well. At the same time, waves are localized well and can be attributed to individual source features.

4 Assessment

In this section we use the methods described in Sect. 3 and follow the assessment approach outlined in Sect. 2 in order to quantify to which accuracy GWMF can be inferred from MLT limb observations and how many independent across-track points (i.e., how many measurement tracks) are required. The assessment measure is the comparison of global GW momentum flux values from simulated temperature observations to the values directly obtained from wind fluctuations of the full model fields. In addition, we consider how well spectral information is conserved when reducing the number of observation tracks. The assessment starts from consistent fluctuations in temperature and wind velocities after scale separation applied to a single snapshot of full model data (see Sect. 2.1 and Fig. 2); i.e., the scale separation is applied but not subject to the assessment.

4.1 Polarization relations

In the proposed mission concept GWMF is inferred from 3-D temperature structures. This requires that (a) polarization relations are also valid in the MLT under nonlinear conditions of many GWs approaching a critical level and that (b) a few-wave decomposition approach such as S3D is an adequate method for determining the 3-D wave vectors of the leading GWs (see Question 1 in Fig. 2). This is tested here using the approach from Sect. 2.3 (see the data flow to the upper yellow diamond of Question 1 in Fig. 2).

In order to verify that the momentum flux based on the S3D-derived wave parameters can be correctly calculated, a comparison with momentum flux estimates directly from wind residuals is carried out. As outlined in Sect. 2.2, we use zonal mean GWMF (without the correction for Coriolis force) as a true reference:

$$\begin{array}{}\text{(9)}& ({\mathit{F}}_x,{\mathit{F}}_y)=\mathit{\rho }\left(\stackrel{\mathrm{‾}}{{\mathit{u}}^{\prime }{\mathit{w}}^{\prime }},\stackrel{\mathrm{‾}}{{\mathit{v}}^{\prime }{\mathit{w}}^{\prime }}\right),\end{array}$$

where F_x and F_y are the zonal and meridional components, respectively. Following the derivation of the vertical flux of horizontal gravity wave pseudomomentum by Ern et al. (2004), we express Eq. (9) in terms of the residual temperature wave amplitude $\widehat{T}$, wavenumbers k, l, and m, and intrinsic frequency $\widehat{\mathit{\omega }}$ as

$$\begin{array}{}\text{(10)}& ({\mathit{F}}_x,{\mathit{F}}_y)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }{\left({\displaystyle \frac{g}{N}}\right)}^{\mathrm{2}}{\left({\displaystyle \frac{\widehat{T}}{T}}\right)}^{\mathrm{2}}{\displaystyle \frac{(\mathit{k},\mathit{l})}{\mathit{m}}}\times {\left(\mathrm{1}-{\displaystyle \frac{f^{\mathrm{2}}}{{\widehat{\mathit{\omega }}}^{\mathrm{2}}}}\right)}^{-\mathrm{1}}.\end{array}$$

In Eq. (10), N is the buoyancy frequency, T is the background temperature, and the factor ${\left(\mathrm{1}-f^{\mathrm{2}}/{\widehat{\mathit{\omega }}}^{\mathrm{2}}\right)}^{-\mathrm{1}}$ converts from pseudomomentum to momentum. Equation (10) is a simplified expression omitting two correction terms which are discussed in the supporting material to Ern et al. (2017) and relevant only to high-frequency GWs not considered in this study. The wave parameters of Eq. (10) are acquired through S3D analysis of temperature residuals of the HIAMCM data.

Figure 9 The zonal mean vertical flux of zonal (red) and meridional (blue) momentum. Dashed lines show GWMF calculated from wind residuals (Eq. 9), while solid lines are calculated using wave vectors and residual temperature amplitudes (Eq. 10). Panels (a) and (b) show the fluxes at 130 km altitude, while panels (c) and (d) show fluxes at 75 km. Wind-based fluxes have been smoothed by averaging over 4 km altitude bins centered at the specified altitude and by averaging over 5^∘ latitude bins.

Figure 9a and b show the zonal and meridional components of GWMF at an altitude of 130 km. The temperature-derived flux, based on wave parameters from global evenly distributed S3D analysis volumes, is in good agreement with the wind-based method used for reference. There is dissimilarity in the zonal GWMF at around 50^∘ N, indicating a region of strong wind shear and/or a breakdown of the validity of the linear approximation. Likewise, at lower altitude the methods are consistent, with some differences in the zonal component around 25^∘ S (Fig. 9c) and in the meridional component across the Southern Hemisphere midlatitudes (Fig. 9d), suggesting dynamics that are not completely captured. As a whole, the results from the two different methods are in good agreement and confirm that Eq. (10), based on wave properties from S3D analysis, is suitable for studying wave dynamics in the MLT.

4.2 How few tracks are required?

One of the major questions concerning our approach is how many measurement tracks are necessary to sufficiently derive the GW parameters in the MLT region within the framework of the proposed two-CubeSat observation strategy (see Question 2 in Fig. 2). For this purpose we test the E2E simulations with a series of a varying number of tracks and swath widths. The detailed results are presented in this section with a focus on the interpretation and intercomparison of GW vectors and zonal mean momentum flux values. All results presented in this subsection are either from sampled orbit data or full E2E simulations (see Sect. 2.3 and the data flow to the middle yellow diamond of Question 2 in Fig. 2).

4.2.1 Analysis of wavelength spectra

Spectra of PGWMF as a function of horizontal and vertical wavelengths are considered for two reasons. First, the choice of the cube size restricts the long wavelength limit of the analyzed spectrum. As described in Sect. 3.7, the desired wavelengths should ideally be in a range of [$\mathrm{1}/\mathrm{2}$, 3] times the cube size, and wavelengths larger than 3 times the cube size are rejected. We hence need to verify that our cube size choice does not cut off major parts of the spectrum. Second, the wavelength spectrum should remain (largely) unchanged when reducing the number of tracks.

Figure 10 Distribution spectra of GW pseudomomentum flux versus logarithmic horizontal and vertical wavenumbers at 75 km (bottom panels), 85 km (middle panels), and 95 km (upper panels) altitude. Varying cube sizes in the across-track direction of 15 tracks (a–c), 5 tracks (d–f), 4 tracks (g–i), and 2 tracks (j–l) are applied. The 15- and 5-track analyses are for sampled data, while 4- and 2-track analyses are for retrieval data with realistic noise. Flagged-out data due to insufficient fitting quality result in more blank bins for 95 km altitude. The gray reference lines through the plots indicate 1000 km horizontal and 10 km vertical wavelength, respectively.

We start with what we deem a good initial value for the cube size: from the HIAMCM set-up the shortest horizontal wavelength is around 156 km, and because from the observational filter of a limb sounder GWs only with horizontal wavelengths longer than 100 km are captured, we expect the shortest horizontal wavelengths of O(200 km). Small-scale GWs with horizontal wavelengths shorter than 200 km are thus not considered here. In order to gain the full picture, satellite observations need to be combined with ground-based systems (e.g., Shiokawa et al., 2009; Nishioka et al., 2013; Chum et al., 2021) observing such shorter scales (see Sect. 5). Based on previous experience and also using a separation scale of zonal wavenumber 18, the longest wavelengths to be considered are O(2000 km). An average vertical wavelength around 12 km was found from SABER data (Ern et al., 2018) for 80 km altitude. Therefore, an initial cube size of 600 km along-track × 420 km across-track × 15 km altitude is selected for fitting the wave vectors and a reduced vertical size of 5 km for refitting the wave amplitudes. With an atm-grid sampling of 30 km × 30 km × 1 km (along-track × across-track × vertical; see Sect. 3.6), this corresponds to 21 × 15 × 15 points. Spectra in terms of horizontal and vertical wavelengths for these initial cube size are shown in Fig. 10a–c for altitudes of 75, 85, and 95 km.

The spectral peak appears at around 600 to 800 km horizontally and 10 to 16 km vertically at all altitudes. All spectra are cut off at longer wavelengths of around 2100 km horizontally and 45 km vertically as the detection upper limit. It results from the limits when filtering reliable fits, which are up to ∼ 3 times the cube size for both horizontal and vertical wavelengths. This, however, does not remove major parts of the spectrum. The spectrum is slightly truncated at shorter horizontal wavelengths around 150 km at an altitude of 75 km, which is not the case for 95 km as waves with longer wavelengths can propagate higher. These wavelength spectra confirm our expectation that, at least for the model data, the target spectral range is covered well by the selected cube size parameters.

In order to investigate the impact of fewer tracks, the number of measurement tracks is reduced in several steps. The proposed mission deploys two SHIs, both with split interferograms to provide four measurement tracks in total. The distance of the track pair from one SHI is assumed to be the distance of the geometric center of the half-interferogram (30 km). The distance between the two satellites can be adjusted by pointing. We assume a gap in the center which can be used to widen the total covered region. In order to study the influence of the gap as well, we use an equidistant sampling of five points here, as well as a four-track configuration in which the center track is not used. If only one SHI instrument were operated, only two measurement tracks would be available. In order to probe all these options, a series of S3D wave analyses were conducted with cube sizes of five, four, and two across-measurement tracks on both sampled and retrieved HIAMCM temperature residual data. The five-track case is included in the simulation and discussed here as it serves as a bridge between the odd and even tracks and offers an opportunity to examine whether the S3D analysis can be performed normally on a reduced number of tracks. The along-track and vertical sizes of the analysis cubes are kept the same, i.e., 21 points and 15 points, respectively. The initial 15-track case is used as a reference for intercomparisons. We have evaluated all cases for both sampled and E2E data, but exemplarily show 15 tracks (Fig. 10a–c) and 5 tracks (Fig. 10d–f) for sampled data and 4 tracks (Fig. 10g–i) and 2 tracks (Fig. 10j–l) for full E2E data.

All three cases of reduced track numbers display the same major patterns as the reference case. It is noted that for the three cases, the spectra also have a horizontal wavelength cutoff limited to around 2100 km in order to provide an upper limit consistent with the reference case. All distributions have their maximum in horizontal wavelength at around 700 km and extend with strong amplitudes to around 300 km. This is compliant with HIAMCM having a Nyquist wavelength of ∼ 100 km and the shortest well-resolved wavelengths of the order of 250 km. For 75 km altitude, there is a secondary peak around 150 km, still inside the Nyquist limit. The fact that the wavelength limit to the short horizontal wavelength side varies with the number of tracks (and hence cube size) is due to an implementation detail of the nested interval variational fit approach. This, however, has no influence on the spectral distribution in the main part of the spectrum. There is a general tendency for larger PGWMF with fewer tracks, which is the highest for the two-track data, which look somewhat blurred, though. In general, all track combinations are suited to recover the spectrum.

4.2.2 Analysis of phase speed and wave direction

A physically more interesting test involves spectra of PGWMF versus ground-based phase speed and direction. The direction is determined from the horizontal wave vector, the intrinsic phase speed can be calculated from the dispersion relation, and the ground-based phase speed is then calculated by Doppler-shifting the phase speed with the large-scale winds. The according PGWMF distribution is shown in Fig. 11 for latitudes of 40 to 70^∘ N. The panels are shown for the same number of tracks and the same altitudes as in Fig. 10. The four cases of different track numbers show common characteristics for both the phase speed values and wave directions. At 75 km altitude there are two lobes, one towards north-northeast (NNE) and a second towards southwest (SW). Maximum phase speeds in these lobes are around 50 m s⁻¹. In addition, there is a widespread background of waves propagating (ground-based) to the east. At 95 km altitude, all low-phase-speed waves are strongly attenuated, and the north-northwest (NNW) lobe is completely removed, presumably by critical-level filtering. The fast eastward waves now prevail. The two-track data seem to have wider spread and some additional features (e.g., fast waves to the southeast, SE) which may be misinterpretations of waves. This is more pronounced at 75 km at the lower edge of the emission layer, where the noise level is higher. In general, however, all cases reproduce the same salient features.

Figure 11 Polar plots of phase speed and direction versus GW pseudomomentum flux in the northern middle and higher latitudes of 40 to 70^∘ at an altitude of 75 km (bottom panels), 85 km (middle panels), and 95 km (upper panels) in line with the four cases as in Fig. 10. The direction is in azimuth angle by which 0^∘ is eastward and 90^∘ is northward. The white dashed radius lines indicate various phase speed values in units of meters per second (m s⁻¹).

Figure 12 Scatter plots of wave direction of the first wave component for 5 tracks (a–c), 4 tracks (d–f), and 2 tracks (g–i) on the y axis versus the reference of 15 tracks on the x axis at an altitude of 75 km (c, f, i), 85 km (b, e, h), and 95 km (a, d, g). A black 1 : 1 identity line is added for the comparison. Note that points around (0, 360^∘) and (360, 0^∘) appear simply due to a mapping of 360^∘ on the y or x axis in plotting.

The deviation of wave direction caused by reducing the number of measurement tracks is further examined by scatter plots of wave directions for the various track-number cases (y axis) against the 15-track reference (x axis) in Fig. 12. The distribution is for the whole globe. We find two clusters of main propagation, a larger one around 180^∘ and a smaller one around 0^∘. For ideal fits, we expect identity; i.e., the same leading waves with the same wave directions would be identified independent of the number of tracks and independent of imposed noise. In this case, all points would be on the 1 : 1 identity line. Indeed, we find that most points cluster around the identity line. There are interesting deviations, though. There are smaller clusters around (0, 180^∘), (180, 0^∘), (360, 180^∘), and (180, 360^∘), which indicate direction flips. The number of direction flips increases with fewer tracks. For the two-track data there is a general loss of ability to determine the propagation direction, which is expressed in vertical stripes around the preferred propagation direction centers. Again, the loss of direction information is most pronounced at 75 km altitude.

4.2.3 Analysis of zonal mean GW momentum flux

Zonal mean GWMF and its vertical gradient constitute the primary goal of the mission and the most stringent test we can apply in the assessment. The dynamical driving and hence the overall structure of the MLT are largely governed by acceleration of the large-scale wind due to GW dissipation. Studies of equatorial oscillations, such as mesospheric semiannual oscillation (MSAO) and mesospheric quasi-biennial oscillation (MQBO), of the general mean circulation and the temperature structure as well as of the formation of an elevated stratopause after a sudden stratospheric warming are largely based on zonal mean GW activity and would highly benefit from accurate PGWMF estimates. This is further explicated is Sect. 5. Furthermore, zonal mean GWMF can be inferred as a true reference from the winds directly as introduced in Sect. 2. We discuss the influence of the observation method on this primary observation aim.

Figure 13 Zonal (a–d) and meridional (f–i) zonal mean GW momentum flux in 1^∘ latitude bins from sampled and retrieved orbit-track HIAMCM temperature data for 1 January 2016 at 06:00 UT using the S3D fitting method and polarization relations. Four cases as in Fig. 10 are listed. GWMF values (e, j) directly inferred from wind fluctuations are also given, which are running-averaged over 5^∘ latitude bins and vertically running-averaged over 5 km altitudes.

Figure 13a–d and f–i depict altitude–latitude cross-sections of the zonal average GWMF for 1 January 2016 at 06:00 UT (i.e., winter in the Northern Hemisphere and summer in the Southern Hemisphere) in 1^∘ latitudinal bins for the four cases of different numbers of observation tracks. Both zonal and meridional components are given. Except for the lowermost row, GWMF is inferred from temperature residual data via S3D wave analysis and polarization relations. We use one snapshot of HIAMCM data but sample with a series of orbits corresponding to 1 week of measurements. In this way we separate sampling issues from general limitations of the method.

Global distribution of GWMF values directly inferred from wind fluctuations as in Eq. (9) is presented in Fig. 13e and j, which serves as a reference. A running average over 5^∘ latitude bins and 5 km altitude bins is applied. All main structures are recovered by the simulated observations.

The main features of Fig. 13 show that GWs intrinsically propagate against the prevailing zonal wind, i.e., eastward propagation in the summer mid-mesosphere around 70 km and westward propagation in the winter hemisphere. The MLT is characterized by strong wind gradients and according wave dissipation and critical-level filtering. This is expressed by strong gradients of GWMF and the fact that the absolute values of momentum flux change by 2 orders of magnitude, and partly, in the zonal direction by a reversal of the direction of zonal GWMF. This zonal mean behavior is consistent with the phase speed spectra of Fig. 11, which show the filtering of the high-GWMF but slow westward GWs between 75 and 95 km altitude and, accordingly, prevailing fast eastward waves at 95 km altitude. This pattern is captured by all four cases of different track numbers.

Figure 14 Line comparison plots of zonal (a–d) and meridional (e–h) zonal mean GW momentum flux calculations in 1^∘ latitude bins from S3D analysis (blue line) and directly from wind fluctuations (orange line) for an altitude of 75 km. It is in line with the four cases as in Fig. 10. GWMF from wind residuals is running-averaged over 5^∘ latitude bins and vertically running-averaged over 5 km altitudes.

Figure 15 Same as for Fig. 14 but for an altitude of 85 km.

Figure 16 Same as for Fig. 14 but for an altitude of 95 km.

For a more quantitative comparison, line plots of momentum flux values from the four cases together with the values computed from wind residual data are presented in Figs. 14, 15, and 16 for altitudes of 75, 85, and 95 km, respectively. The line plots indicate good agreement of the momentum flux values inferred from the observed temperature residuals with the reference from model wind fluctuations. Some deviations are found in the midlatitudes for zonal GWMF components, especially at an altitude of 75 km, at 30^∘ S and 40^∘ N. For meridional GWMF components, discrepancies mainly appear in the Southern Hemisphere at 40–60^∘ S for altitudes of 85 and 95 km. These differences are due to strong vertical gradients and could either be a problem of the wave fitting method to identify the local vertical wavelength associated with the cube center altitude or an effect of nonlinearity and limitations to the linear GW physics employed to calculate GWMF.

In conclusion, the assessment results show that wave analysis of tomographic temperature observations with few observation tracks (down to four tracks and even two tracks) is suitable to gain reliable zonal means of zonal and meridional GWMF.

4.3 How much noise can we afford?

The E2E assessment performed in Sect. 4.2 indicates the viability of the proposed mission concept based on our best estimate of instrument performance. In order to investigate whether further miniaturization of the instrument and related decrease in the signal-to-noise ratio would be feasible (see Question 3 in Fig. 2), the best-estimate noise level superposed on the synthetic spectra is scaled by multiples of 2, 4, 6, 16, and 32. After that, temperature retrieval and S3D wave analysis are performed as before. We focus on the four-track data and assess the E2E results by zonal mean zonal GWMF values (see Sect. 2.3 and the data flow to the lower yellow diamond of Question 3 in Fig. 2).

Figure 17 Line plots of the zonal components of zonal mean GW momentum flux in 5^∘ latitude bins from S3D analysis on retrieved temperature residuals with estimation of different noise levels of 1 (orange line), 2 (blue line), 4 (green line), 8 (cyan line), 16 (magenta line), and 32 (yellow line) at an altitude of 75 km (a), 85 km (b), and 95 km (c). A four-track case is applied for this analysis.

We show in Fig. 17 color-coded line plots for noise levels of 1, 2, 4, 8, 16, and 32 for altitudes of 75, 85, and 95 km, respectively. The resulting distributions displayed in Fig. 17 are noisier and coarser since only a single day of orbits was used. The zonal mean GWMF plots indicate that the GW structure tends to be damped with increasing noise level, but the main features are retained up to a noise level of 8 times the original. For a 32-fold noise level, the wave signals are no longer discernible. Further detailed analysis of zonal mean GWMF for a 4-fold noise level is given in Appendix C. In general, our approach could tolerate enhanced noise up to a factor of 4 higher than our best estimate of the instrument of Kaufmann et al. (2018), which allows for further downsizing of the instrument or higher detector temperatures resulting in higher dark current levels and shot noise.

5 Discussion of science applications

Which scientific questions could be directly addressed if an instrument such that as discussed in this paper were in orbit? And for which studies would we need further ancillary data? A two-step processing chain as it would be applied to in-orbit data is sketched in Figs. 18 and 19. The first part relies only on data directly obtained by retrieval and wave analysis, and the second would involve atmospheric winds that can be calculated, e.g., by data assimilation.

Figure 18 Processing chain including retrieval and GW analysis. Up to a climatology of GWMF in terms of, e.g., zonal means or monthly mean 3-D global distributions (Level 3 data) only data observed by the instrument are required.

Figure 19 Processing chain for scientific interpretation of propagation, source identification, and critical-level filtering. In order to determine the ground-based frequency and perform ray-tracing background wind velocities for the observation location and time are required. These can be generated as Level 4 data via assimilating various data sets (including the observations described here) in a consistent manner.

Gravity wave momentum flux can be calculated using only the observations made by the proposed instrument (Fig. 18). From the GWMF values of single GW events, climatological distributions, such as average monthly mean maps, zonal means, drag from vertical gradients of GWMF, and spectral distributions in terms of horizontal and vertical wavelength or intrinsic phase speed, can be generated. Such distributions are always the starting point of more in-depth scientific investigations and can be used for interpretations of potential sources. Together with climatologies of winds and temperatures from other observations covering the MLT altitude range (e.g., the URAP climatology; Swinbank and Ortland, 2003) they can be used to gain a first assessment of the momentum balance. Finally, spectral distributions including the direction can be used to distinguish different pathways into the MLT (as outlined in the introduction) and hence provide the information to distinguish between conceptually different modeling approaches (both GW-allowing GCMs, which still contain tunable parameters, and explicit physical GW models) better than the GW variances and absolute values of GWMF we can use nowadays.

Even more investigations are facilitated if estimates of the large-scale winds at observation time are also available. In the stratosphere such wind data are regularly generated by assimilation systems operated by numerical weather prediction centers. In the mesosphere geostrophic winds are often used, and zonal mean values were generated up to 90 km altitude (Ern et al., 2013b; Smith et al., 2017). In order to gain a 3-D and time-dependent picture including tides, data assimilation is a powerful tool (e.g., Eckermann et al., 2009; Pedatella et al., 2018, 2020).

Zonal mean winds already allow assessment of the driving of large-scale wind patterns by GWs. Examples for these are studies based on absolute values of GWMF from limb scanning. These are the most reliable estimate of global GWMF distributions we can currently gain, but they are limited because of the lack of direction information. In zones of strong vertical wind shear and an environment wherein slow-phase-speed GWs dominate, we can assume that the vertical gradient of the momentum flux corresponds to drag directed opposite to the shear (i.e., by GWs causing the shear layer to propagate downward). This concept has been very successfully used, e.g., in studies of the stratospheric quasi-biennial oscillation (Ern et al., 2014), the excitation of quasi-2 d waves (Ern et al., 2013b), and the role of GWs in sudden stratospheric warmings (Ern et al., 2016). However, in regions where substantial filtering of slower-phase-speed waves has already occurred at lower altitudes and only faster phase speeds survive or when GWs dissipate primarily due to the increasing amplitudes the waves attain when they propagate upwards into regions of lower density, such simplifications cease to work. That can be seen, for instance, in the discussion of the mesospheric semiannual oscillation (Ern et al., 2015, 2021) wherein arguments for the direction of drag became complex and indirect. This increased complexity is the mark of many of the large-scale global wind patterns in the mesosphere.

New 3-D data are a great step forward from existing observations: in the MLT region, currently existing satellite instruments do not provide 3-D information about observed GWs such that directional GWMF and directional GW drag cannot be derived. Conventional limb sounders, such as SABER or the Microwave Limb Sounder (MLS), provide only a single measurement track of altitude profiles. Accordingly, these observations are limited to GW variances (e.g., Jiang et al., 2005; Hocke et al., 2016) and GW absolute momentum fluxes (e.g., Ern et al., 2018, 2022). Solar occultations, for example by the Solar Occultation for Ice Experiment (SOFIE), are even more limited due to their sparse sampling and limited global coverage (e.g., X. Liu et al., 2014; Thurairajah et al., 2014). Other satellite instruments provide 2-D horizontal information, but GW vertical wavelengths cannot be determined (e.g., Rong et al., 2018; England et al., 2020). This means that a climatology of directional GW momentum fluxes is currently missing in the MLT region.

The proposed instrument combined with background wind information would resolve ambiguities such as those encountered by Ern et al. (2015, 2021) immediately by combining the observed GW parameters with estimates of background winds (see Fig. 19). For this, a relatively coarse spatial resolution including only wavenumbers 0–6 and a moderate altitude resolution would be sufficient. For instance, calculating the Doppler shift allows converting from intrinsic to ground-based phase speed diagrams. The ground-based phase speed spectrum could then directly be compared to the wind velocities by means of blocking diagrams as introduced by Taylor et al. (1993) and the type of GW dissipation (critical level versus saturation of growing amplitudes) inferred. The vertical gradient of the zonal or regional mean of GWMF can be used to infer net drag including its direction and taking into account GWs from all propagation directions. Also, in complex situations such as the mesospheric SAO, the contribution of GWs to the driving can be directly estimated. For such applications zonal mean winds or coarse spatial representations would suffice and a general guidance of the GCM by the assimilation system would produce novel insights.

It should be noted here that current-day assimilations still have difficulties in particularly capturing those dynamical features wherein GWs play a major role (Harvey et al., 2022a). However, when more sensors gain direction-resolved GW information, a stronger need for such fields hopefully drives improvements. The challenge in the MLT is the superposition of large-scale flow and tides, which can have similar wind amplitudes as the magnitude of the large-scale flow. However, using geostrophic winds, tidal determination from observations, and dedicated tidal models (e.g., Nguyen and Palo, 2013; Pedatella et al., 2016) there are alternatives which should also suffice for most of the investigations sketched here.

A more sophisticated way to study the interaction of GWs and large-scale winds is forward and backward ray tracing. This is more accurate than blocking diagrams as, for instance, lateral propagation and avoidance of critical levels are taken into account (Preusse et al., 2009a; Kalisch et al., 2014; Thurairajah et al., 2014; Chen et al., 2019; Thurairajah et al., 2020). Furthermore, tides, e.g., diurnal and semi-diurnal tides, cause changes in the large-scale winds at similar timescales as the periods of the GWs propagating through these winds. This causes refraction of the ground-based wave frequency (Senf and Achatz, 2011), thus also allowing waves, which would be expected to be filtered by a tide based on the original ground-based frequency at the source level, to propagate further. Having identified critical levels in backward ray tracing, however, this is a clear sign for in situ generation of GWs in the middle atmosphere by unstable jets or secondary wave generation. Also, for such studies a coarser representation of the background winds is sufficient.

Most demanding is the backward tracing of observed waves to potential tropospheric sources. Backward or forward ray tracing of individual wave events requires a full wave characterization and for temperature observations hence requires three dimensions. Backward ray tracing from 3-D data was used in previous studies for source identification. Examples for mesoscale stratospheric waves are studies by Preusse et al. (2014), Krisch et al. (2017), Perrett et al. (2021), Strube et al. (2021), and Geldenhuys et al. (2021). High-frequency waves in the MLT have also been studied by such means (Wrasse et al., 2006; Pramitha et al., 2015). Backward ray tracing is a comprehensive method to infer whether GWs can be followed down to the troposphere or meet a nontransparent level higher up – the latter indicates a source in the middle atmosphere, for instance by secondary wave generation. This is an advance over simply considering phase speeds at the observation and wind fields below as all the spatial and temporal changes in the wave parameters along the ray path are taken into account. Whether the precision of the wave parameters and the accuracy of the background winds from assimilation will be sufficient to trace mesoscale waves back to individual tropospheric sources will need to be studied in more detail. Upward ray tracing can then be used to investigate the interaction of GWs with the background and to provide a complementary drag assessment to the vertical gradient of PGWMF. For instance, which waves can reach the MLT could be assessed.

Although it is not sure whether global wind observations, or reliable winds from data assimilation, will be available in the upper mesosphere–lower thermosphere at the time the instrument will be in operation, the gravity wave data set that we expect to obtain from this novel observation method will be quite unique and of great value in itself. Even without wind observations, studies based on the observed directional gravity wave momentum flux can be performed in a climatological sense, for example by comparison with zonal wind climatologies or climatologies of atmospheric tides. Particularly, the interaction between gravity waves and tides is not well understood and offers a wide field of applications. Further, the novel gravity wave data set can be used to identify cases of excitation of secondary gravity waves. This can be performed by identifying fishbone-like structures in along-track and altitude cross-sections without the need for background wind information (e.g., Vadas et al., 2018; Vadas and Becker, 2018). For these kind of studies, the relatively short along-track sampling of 30 km, combined with a tomographic retrieval, will be very beneficial.

In the context of momentum balance studies, the observational filter of a limb sounder needs to be taken into account. The observational filter allows observing only waves with horizontal wavelengths longer than O(100 km). It is evident that this misses an important part of the GW spectrum. For instance, ground-based systems (e.g., Shiokawa et al., 2009; Nishioka et al., 2013; Chum et al., 2021) observing such shorter scales indicate large GWMF values which, upscaled to global distributions, would indicate that a larger part of the GWMF is contributed by the shorter scales. However, ground stations are situated only on land, often in hot-spot regions, rather than over open oceans, which biases the distributions. High-altitude GCMs without GW parameterization provide an argument for the importance of longer scales. They are even more restricted to long horizontal wavelengths than limb-sounding observations (e.g., Sato et al., 1999, 2012; Siskind, 2014; Becker and Vadas, 2020). It is one of the puzzles of atmospheric dynamics that these GCMs still produce a realistic atmosphere even above the stratopause where GWs become increasingly important. Further indications of the importance of scales visible to a limb sounder are the momentum budget studies discussed in the Introduction. Still, the importance of different scale GWs for the driving of the upper mesoscale is an unsolved puzzle.

For the stratosphere nowadays global simulations with a grid distance of a few kilometers and even down to 1 km are available (e.g., Stephan et al., 2019a, b; Polichtchouk et al., 2022), which allow assessment of how large a fraction of GWMF is likely to be missed due to the observational filter. However, such simulations have an upper altitude limit of less than 80 km and most of them damp GWs in a strong sponge layer above 40 km (Preusse et al., 2014, and references therein). The best estimate for such lower altitudes is that GWMF is about equally partitioned to scales shorter and longer than 100 km horizontal wavelength, respectively. The almost vanishing influence of the observational filter in our study is thus due to the fact that current state-of-the-art models do not resolve the shorter scales, and it is reasonable to assume that in reality only roughly half of the GWMF would also be observed in the mesosphere. The partitioning of momentum flux between different wavelength regions is, however, one of the major unknowns in the field of GWs, and low-error, well-characterized global observations of GWMF would be an essential and necessary step forward to answer this fundamental question. A brief outline of the current state of knowledge is given in Appendix D. In the lower thermosphere short-period, short-horizontal-wavelength GWs become strongly damped (e.g., Pitteway and Hines, 1963; Vadas, 2007; Yigit et al., 2009), and hence likely an even larger fraction of that part of the GW spectrum for which GWs can propagate upward is quantified.

6 Conclusions

The energy and momentum gravity waves deposit in the mesosphere and lower thermosphere are major drivers of the dynamics and hence the entire structure of this region. In order to understand these processes, we require global observations of gravity wave pseudomomentum flux (PGWMF), its vertical gradient, and its spectral composition. In particular, the distribution of PGWMF with respect to phase speed and direction allows direct inference of the interaction of GWs with the background winds and the mechanisms of how they drive global-scale patterns such as the general mean circulation and tropical oscillations. By observing 3-D temperature structures, one should be able to gain such information from space. The current study investigates aspects of the concrete realizations and whether it is sufficient to sample gravity waves in that part of the atmosphere with only two to four instantaneous fields of view.

To assess this question, an end-to-end simulation comprising orbit simulation, temperature retrieval, and 3-D wave analysis was performed, and the performance was assessed on the basis of GWMF comparisons. We considered zonal mean values of zonal and meridional GWMF deduced from the simulated observations and compared this to a reference distribution deduced directly from the model winds. In addition, we investigated phase speed spectra and the ability to characterize the direction of single wave events.

Our simulations show that zonal mean values are most robust and that the direction of wave events is most sensitive with respect to both noise and a low number of observation tracks. In general, a four-track instrument with nominal performance provides reliable wave quantification and can be employed for both the generation of zonal mean climatologies and studies of the interaction of GWs with the background winds.

For reasons of computational efficiency, we made the following assumptions: since the HIAMCM resolves only GWs longer than about 200 km horizontal wavelength, we have to rely on a horizontal background removal based on global-scale wave estimates. Accordingly, scale separation is not included in the assessment process. In addition, we assume a constant track distance of the split interferogram of 30 km (which implicitly assumes a constantly lit scene). Such assumptions are necessary to make the study feasible and are, on the whole, justified, but they tend to overestimate the reliability of the results somewhat. An important finding is that the signal-to-noise levels of the temperature data assumed in this study can be relaxed by a factor of 4 to allow for further miniaturization of the observation system. For the gravity wave spectrum used in this study, two fields of view parallel to the orbital track are sufficient to resolve most of the spectral distribution, whereas the direction of individual GWs cannot be resolved as well as in simulations assuming four or more tracks.

Summarizing, the most important conclusions obtained are listed below.

• Assumptions made for the S3D fitting technique are non-critical in practical application.

• Polarization relations are valid up to 130 km altitude.

• Four closely spaced fields of view are sufficient to resolve the spectral distribution and the direction of gravity waves in the mesosphere. Two fields of view are sufficient to resolve the spectral distribution, but uncertainties of inferred GW direction information are increased.

• Inferred distributions are still meaningful even if noise is increased by a factor of 4 with respect to the assumed instrument performance.

Appendix A: Photochemical reaction parameters used in the O₂ airglow emission simulation

Following Sheese (2009), the A-band volume emission rate is given by

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}\mathit{\eta }& =(P_{\text{A band}}+P_{\text{B band}}+P_{{\mathrm{Phot}}_{{\mathrm{O}}_{\mathrm{3}}}}+P_{{\mathrm{Phot}}_{{\mathrm{O}}_{\mathrm{2}}}}\\ & +P_{\mathrm{Barth}}){\displaystyle \frac{A_{\mathrm{762}}}{P_{\mathrm{loss}}}},\end{array}\end{array}$$

where the first term describes the production rate of the excited state ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$, and the second term is the fraction of the O₂ A band compared to all losses of the excited state ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$. A₇₆₂ is the Einstein coefficient of the A band and can be calculated by A₇₆₂=F_cA_1Σ, where F_c is the A-band Franck–Condon factor, and A_1Σ is the Einstein coefficient of all radiative transitions from the excited state ${\mathrm{O}}_{\mathrm{2}}(b^{\mathrm{1}}{\mathrm{\Sigma }}_g^{+},v=\mathrm{0})$ corresponding to Reaction (R19) in Table A1. In the following, [⋅] corresponds to the number densities of molecules. The loss component can be calculated by

$$\begin{array}{}\text{(A2)}& P_{\mathrm{loss}}=A_{\mathrm{1}\mathrm{\Sigma }}+k_{\mathrm{0}}\left[{\mathrm{N}}_{\mathrm{2}}\right]+k_{\mathrm{3}}\left[{\mathrm{O}}_{\mathrm{3}}\right]+k_{\mathrm{4}}\left[{\mathrm{O}}_{\mathrm{2}}\right]+k_{\mathrm{6}}\left[\mathrm{O}\right],\end{array}$$

where k₀, k₃, k₄, and k₆ correspond to the quenching rates of Reactions (R20)–(R23) in Table A1. The A-band photo-absorption is calculated by

$$\begin{array}{}\text{(A3)}& P_{\text{A band}}=g_{\mathrm{A}}\left[{\mathrm{O}}_{\mathrm{2}}\right],\end{array}$$

where g_A is the photochemical reaction coefficient of Reaction (R1) in Table A1. The B-band photo-absorption is calculated by

$$\begin{array}{}\text{(A4)}& {\displaystyle }{P}_{\mathrm{B}\mathrm{band}}={\displaystyle \frac{K{g}_{\mathrm{A}}\left[{\mathrm{O}}_{\mathrm{2}}\right]}{A_{\mathrm{771}}+K+K_{\mathrm{3}\mathrm{B}}\left[{\mathrm{O}}_{\mathrm{3}}\right]}},\text{(A5)}& {\displaystyle }K=K_{\mathrm{0}B}\left[O\right]+K_{\mathrm{1}B}\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_{\mathrm{2}B}\left[{\mathrm{N}}_{\mathrm{2}}\right],\end{array}$$

where K_0B, K_1B, K_2B, and K_3B are the quenching rates of Reactions (R4)–(R7), g_B is the photochemical reaction coefficient of Reaction (R2), and A₇₇₁ is the Einstein coefficient of Reaction (R3) in Table A1.

Bucholtz et al. (1986)Bucholtz et al. (1986)Rothman et al. (2013)Yankovsky and Manuilova (2006)Yankovsky and Manuilova (2006)Yankovsky and Manuilova (2006)Yankovsky and Manuilova (2006)Sheese (2009)Sheese (2009)Sheese (2009)Rothman et al. (2013)Green et al. (2000)Sander et al. (2011)Sander et al. (2011)McDade et al. (1986)McDade et al. (1986)McDade et al. (1986)Nicholls (1965)Rothman et al. (2013)Sander et al. (2011)Sander et al. (2011)Sander et al. (2011)Sander et al. (2011)

Table A1 Photochemical reaction parameters used in the simulation of the A-band production and loss mechanisms.

The production rate due to the photolysis of O₂ and O₃ can be calculated by

$$\begin{array}{}\text{(A6)}& {\displaystyle }{P}_{{\mathrm{Phot}}_{{\mathrm{O}}_{\mathrm{2}}}}={\displaystyle \frac{(J_{\mathrm{SCR}}+J_{{\mathrm{L}}_{\mathit{\alpha }}})\left[{\mathrm{O}}_{\mathrm{2}}\right]\mathit{\phi }{k}_{\mathrm{1}}\left[{\mathrm{O}}_{\mathrm{2}}\right]}{A_{\mathrm{1}D}+k_{\mathrm{1}}\left[{\mathrm{O}}_{\mathrm{2}}\right]+k_{\mathrm{2}}\left[{\mathrm{N}}_{\mathrm{2}}\right]}}\text{(A7)}& {\displaystyle }{P}_{{\mathrm{Phot}}_{{\mathrm{O}}_{\mathrm{3}}}}={\displaystyle \frac{J_{\mathrm{H}}\left[{\mathrm{O}}_{\mathrm{3}}\right]\mathit{\phi }{k}_{\mathrm{1}}\left[{\mathrm{O}}_{\mathrm{2}}\right]}{A_{\mathrm{1}D}+k_{\mathrm{1}}\left[{\mathrm{O}}_{\mathrm{2}}\right]+k_{\mathrm{2}}\left[{\mathrm{N}}_{\mathrm{2}}\right]}}\end{array}$$

where J_SCR, $J_{{\mathrm{L}}_{\mathit{\alpha }}}$, and J_H are the photolysis coefficients of Reactions (R8), (R9), and (R10), A_1D is the Einstein coefficient of Reaction (R11), and k₁ and k₂ are the quenching rates of Reactions (R13) and (R14) in Table A1. J_H is given by Sheese (2009). J_SCR and $J_{{\mathrm{L}}_{\mathit{\alpha }}}$ are calculated by following Sheese (2009). Hereby, the Lyman-α absorption cross-section and absorption quantum yield are set to 1.0 × 10⁻²⁰ cm² and 0.55, respectively (Reddmann and Uhl, 2003). Regarding the Schumann–Runge continuum, the absorption cross-section is given by Yoshino et al. (2005) over the given wavelength range and the absorption quantum yield is set to 1 (Brasseur and Solomon, 2005). The solar flux at the top of the atmosphere is taken from the SORCE Solar Spectral Irradiance (SSI) data product and can be obtained from https://lasp.colorado.edu/home/sorce/data/ (last access: 5 December 2022).

The production rate due to the Barth process can be calculated by

$$\begin{array}{}\text{(A8)}& P_{\mathrm{Barth}}={\displaystyle \frac{k_{\mathrm{5}}\left[\mathrm{O}{]}^{\mathrm{2}}\right[M\left]\right[{\mathrm{O}}_{\mathrm{2}}]}{C_{{\mathrm{O}}_{\mathrm{2}}}\left[{\mathrm{O}}_{\mathrm{2}}\right]+C_{\mathrm{O}}\left[\mathrm{O}\right]}},\end{array}$$

where k₅ is the quenching rate of Reaction (R15) in Table A1, and $C_{{\mathrm{O}}_{\mathrm{2}}}$ and C_O are the fitting parameters of the simplified model proposed by McDade et al. (1986). $C_{{\mathrm{O}}_{\mathrm{2}}}$ describes the relative quenching rates of ${\mathrm{O}}_{\mathrm{2}}^{\ast }$ with O₂ corresponding to Reaction (R16), and C_O describes the relative quenching rates of ${\mathrm{O}}_{\mathrm{2}}^{\ast }$ with O corresponding to Reaction (R17) in Table A1, which is a three-step reaction as depicted in Fig. 4a.

Appendix B: Known limitations of the O₂ A-band forward model

To allow for repeated tomographic retrievals of a week of data, the employed forward model makes some simplifying assumptions that need to be revisited before applying it to actual measurements. First, for daytime measurements, the model requires the amount of in-scattered solar light at all positions. This is currently being tabulated based on climatological conditions and not based on actual volume-mixing ratios. Second, the model does not take the finite extent of the line shape into account. The support of the line shape is very small compared to the spectral resolution of the instrument such that this becomes problematic only for the treatment of self-absorption: we simplify this by computing the amount of self-absorption from the center of the line, which is a strict overestimation of actual self-absorption that decreases in the line wings. In this fashion, the access to radiation from lower altitudes (below 80 km) becomes increasingly limited, which we deem a worst-case assumption for the purpose of this study; this is asserted by runs without any self-absorption that give better-quality results down to 60 km altitude. Both issues are straightforward to address by performing the required computations.

Appendix C: Zonal mean distribution of GWMF for a 4-fold increased noise level

Illustrated in Fig. C1 is the global cross-section distribution of zonal mean GWMF in 1^∘ latitude bins, similar to Fig. 13c of the four-track case, but for a noise level enhanced by a factor of 4. The general GW features in Fig. C1 are comparable to Fig. 13c and h, with the wind reversal and strong vertical gradients still clearly visible. In the altitude region above 80 km, these GW structures are in general well preserved in the averaged GWMF distribution, except for some outliers seen in the meridional components at 20–30^∘ N for an altitude of around 85 km. Below 80 km, the latitude bins appear to be biased by noise perturbations. Particularly at 0–30^∘ N and 60–40^∘ S no consistent GW patterns can be identified. It reveals that the wave analysis method becomes more sensitive to noises at lower altitudes. Overall, this GWMF distribution demonstrates that our approach of wave analysis is fairly stable despite increasing the noise level by a factor of 4.

Figure C1 Zonal (a) and meridional (b) components of zonal mean GW momentum flux in 1^∘ latitude bins from retrieved temperature residuals with an amplified 4-fold noise level applied to a four-track case using the S3D fitting method and polarization relations.

Appendix D: The relative importance of different horizontal scales

Coupling of the different layers of the atmosphere is achieved by freely propagating (i.e., internal) GWs; other wave modes such as evanescent waves or waves in a wave guide do not contribute. This limits the horizontal range of GWs which have this potential to horizontal wavelength λ_h in the range of 20–30 km on the short wavelength side and to roughly 1000–2000 km on the long wavelength side (Preusse et al., 2008). Shorter wavelengths exist as well (Fritts et al., 2017), but they are important in the GW cascade when it comes to transferring momentum and energy finally to turbulence, i.e., as the last step of the cascade. A very rough first approximation of the partitioning of GWMF between different scales can be gained by assuming a universal scaling law ${T^{\prime }}^{\mathrm{2}}\propto k^{-\mathrm{5}/\mathrm{3}}$, where the wind or temperature fluctuations scale exponentially with the horizontal wavenumber. This lets us expect that a horizontal wavelength of 100 km roughly partitions GWMF into equal parts. A detailed review by Preusse et al. (2008) found that the available evidence at the time was roughly compliant with such a partitioning. Note that here we discuss the global scale: over localized sources shorter horizontal wavelengths sometimes have larger GWMF by factors (e.g., Kruse et al., 2022), but this is then also connected to a limited area. Which evidence did we gain since the previous review of Preusse et al. (2008)? First, there are now a number of studies using the vertical derivative of absolute GWMF inferred from limb-sounding satellite instruments. It is argued in these studies that in regions of strong wind shear in particular, those GWs propagating opposite to the winds in the lower altitudes break. Despite the fact that the absolute GWMF should be underestimated by the observation technique, the contribution to the total GW acceleration expected from momentum closure is very often more than 50 % (Ern et al., 2013b, 2014, 2015, 2016). Second, we now have more model studies with GW-permitting resolution not using any GW parameterization (e.g., H. L. Liu et al., 2014; Siskind, 2014; Becker and Vadas, 2020; Okui et al., 2021; Becker et al., 2022; Okui et al., 2022). These models reach at least similar performance for the MLT and better performance for the thermosphere than GW-parameterizing GCMs as they include middle-atmosphere sources and oblique propagation, despite the fact that they resolve only GWs with horizontal wavelengths longer than ≈ 200 km. Likely, these models then somewhat overestimate the scales they resolve. Third, there are now GCM runs from numerical weather prediction (NWP) models with a much higher resolution (Stephan et al., 2019a; Polichtchouk et al., 2022). These data also indicate a roughly equal partitioning of GWMF between scales smaller and larger than a horizontal wavelength of λ_h = 100 km. Thus, assuming that a horizontal wavelength larger than 100 km conveys roughly half of the momentum flux is currently the best estimate we have. Unless we perform high-accuracy global GWMF observations this fundamental question will remain at least partly subject to speculation.