the Creative Commons Attribution 4.0 License.
Special issue: The SPARC Reanalysis Intercomparison Project (SRIP) (ACP/ESSD...
Research article 10 Aug 2020
Research article  10 Aug 2020
Lagrangian gravity wave spectra in the lower stratosphere of current (re)analyses
 ^{1}Laboratoire de Météorologie Dynamique (LMD/IPSL), École polytechnique, Institut polytechnique de Paris, Sorbonne Université, École normale supérieure, PSL Research University, CNRS, Paris, France
 ^{a}formerly at: Forschungszentrum Jülich (IEK7: Stratosphere), Jülich, Germany
 ^{1}Laboratoire de Météorologie Dynamique (LMD/IPSL), École polytechnique, Institut polytechnique de Paris, Sorbonne Université, École normale supérieure, PSL Research University, CNRS, Paris, France
 ^{a}formerly at: Forschungszentrum Jülich (IEK7: Stratosphere), Jülich, Germany
Correspondence: Aurélien Podglajen (aurelien.podglajen@lmd.polytechnique.fr)
Hide author detailsCorrespondence: Aurélien Podglajen (aurelien.podglajen@lmd.polytechnique.fr)
Due to their increasing spatial resolution, numerical weather prediction (NWP) models and the associated analyses resolve a growing fraction of the gravity wave (GW) spectrum. However, it is unclear how well this “resolved” part of the spectrum truly compares to the actual atmospheric variability. In particular, the Lagrangian variability, relevant, for example, to atmospheric dispersion and to microphysical modeling in the upper troposphere–lower stratosphere (UTLS), has not yet been documented in recent products.
To address this shortcoming, this paper presents an assessment of the GW spectrum as a function of the intrinsic (air parcel following) frequency in recent (re)analyses (ERAInterim, ERA5, the ECMWF operational analysis and MERRA2). Longduration, quasiLagrangian balloon observations in the equatorial and Antarctic lower stratosphere are used as a reference for the atmospheric spectrum and are compared to synthetic balloon observations along trajectories calculated using the wind and temperature fields of the reanalyses. Overall, the reanalyses represent realistic features of the spectrum, notably the spectral gap between planetary and gravity waves and a peak in horizontal kinetic energy associated with inertial waves near the Coriolis frequency f in the polar region. In the tropics, they represent the slope of the spectrum at low frequency. However, the variability is generally underestimated even in the lowfrequency portion of the spectrum. In particular, the nearinertial peak, although present in the reanalyses, has a reduced magnitude compared to balloon observations. We compare the observed and modeled variabilities of temperature, zonal momentum flux and vertical wind speed, which are related to low, mid and highfrequency waves, respectively. The probability density function (PDF) distributions have similar shapes but show increasing disagreement with increasing intrinsic frequency. Since at those altitudes they are mainly caused by gravity waves, we also compare the geographic distribution of vertical wind fluctuations in the different products, which emphasizes the increase of both GW variance and intermittency with horizontal resolution. Finally, we quantify the fraction of resolved variability and its dependency on model resolution for the different variables. In all (re)analysis products, a significant part of the variability is still missing, especially at high frequencies, and should hence be parameterized. Among the two polar balloon datasets used, one was broadcast on the Global Telecommunication System for assimilation in NWP models, while the other consists of independent observations (unassimilated in the reanalyses). Comparing the Lagrangian spectra between the two campaigns shows that the (re)analyses are largely influenced by balloon data assimilation, which especially enhances the variance at low GW frequency.
Atmospheric gravity waves (GWs) are mesoscale motions with largescale impacts notably through three mechanisms. First, they transport momentum from lower levels and deposit it higher up in the atmosphere, which forces largescale circulations (Andrews et al., 1987), such as the quasibiennial oscillation (QBO; Baldwin et al., 2001). Second, they generate smallscale turbulence (e.g., when breaking), which contributes to mixing atmospheric trace constituents (Podglajen et al., 2017) and diabatic heating. Third, GWs induce temperature and wind fluctuations which impact the formation and microphysical properties of clouds (e.g., cirrus clouds; Potter and Holton, 1995) and aerosols.
Because of those largescale effects, GWs need to be represented in global atmospheric models. In current climate models with resolutions in the order of 100 km, GWs are mostly unresolved and need to be parameterized. In contrast, global weather forecast models, which currently have resolutions down to about 10 km or less, may now start to resolve a significant portion of the GW spectrum (e.g., Preusse et al., 2014; Jewtoukoff et al., 2015; Holt et al., 2016). However, the exact fraction of resolved GWs depends not only on the nominal resolution of the model but also on the parameterized diffusion and on the representation of wave sources like tropospheric convection (Stephan et al., 2019). A common flaw of resolved GWs in models appears to be an underestimation of wave amplitude and an overestimation of horizontal wavelengths (e.g., Geller et al., 2013; Holt et al., 2017). Even in models with realistic GW generation, a lack of realism in the propagation and dissipation of the waves often renders additional GW parameterization necessary in order to obtain a realistic general circulation (Holt et al., 2016, 2017). Hence, for a given model with a given resolution, it is not clear a priori which GWs are represented and which should be parameterized. This issue will become increasingly problematic as the models refine their resolution in the socalled “gray zone”, resolving a larger part (yet not all) of the GW spectrum and its sources.
Crucial for climate and weather forecast models, the question of the fraction of resolved GWs is also important in (re)analyses. Those products have been used to investigate some properties of the GW field (e.g., Preusse et al., 2014). Reanalyses are also widely employed as input to trajectory calculations, notably (but not only) in the upper troposphere–lower stratosphere (UTLS), in order to understand, for example, transport (e.g., Tzella and Legras, 2011), chemistry (e.g., Konopka et al., 2010) and cirrus cloud formation (e.g., Jensen and Pfister, 2004; Schoeberl et al., 2015; Ueyama et al., 2015). Among those three processes, cloud formation (Spichtinger and Krämer, 2013; Dinh et al., 2016) is especially sensitive to mesoscale fluctuations of wind and temperature; thus a number of parameterizations have been developed to account for them (e.g., Bacmeister et al., 1999; Jensen and Pfister, 2004; Gary, 2006). A difficulty is that the parameterized fluctuations need to be similar to the perturbations experienced by air parcels (i.e., Lagrangian).
Recently, Podglajen et al. (2016b) used longduration superpressure balloon (SPB) observations to characterize Lagrangian wind and temperature fluctuations due to GWs in the lower stratosphere. SPBs are especially suited for gravity wave studies since they follow the wind and directly provide access to the intrinsic frequency $\widehat{\mathit{\omega}}$ of the waves. Podglajen et al. (2016b) proposed a parameterization of the vertical wind and temperature fluctuations for Lagrangian trajectory models that use (re)analyses to compute the trajectory path. A version of the parameterization approach that is simpler to implement and its salient features are presented in Kärcher and Podglajen (2019). However, both Podglajen et al. (2016b) and Kärcher and Podglajen (2019), like other authors, implicitly assumed that most of the GWs seen in the observations were absent from the original reanalysis products. In the present paper, we compare GWinduced fluctuations in modern (re)analyses from the same point of view as SPBs, i.e., a Lagrangian point of view, to determine which part of the GW spectrum they resolve and how this resolved part compares with observations. Particular focus will be spent on the intermittency of the gravity wave fluctuations in the (re)analyses.
Besides the interest for studies using Lagrangian trajectories, the comparison presented here also serves another purpose but this time from the point of view of model developers. For different models there are different criteria of success, a priori emphasizing the representation of the largescale circulation. As the resolution of the models increases, a richer array of phenomena is present in the resolved fields, and the question arises of whether or not to assimilate information on these processes. In that respect, assessing model errors in the GW field, which is described in both observations (e.g., radiosondes, GPS temperature profiles) and model output, is essential from the modelers' point of view. The longduration balloon dataset provides a unique opportunity to perform such a task and assess the realism of modeled GWs in the lower stratosphere.
The paper is organized as follows. Section 2 presents the balloon dataset and reanalyses used, as well as the comparison methodology. Then, the results and comparisons are presented in Sect. 3 and discussed in Sect. 4. Finally, summary and conclusions are provided in Sect. 5.
2.1 Longduration balloon observations
Most observations of the atmosphere are obtained either at a fixed location (radar, lidar) or in motion relative to the flow (satellite, aircraft). In both cases, the relative speed between the measuring instruments and the air hampers direct Lagrangian analysis of flow variability. One of the platforms which overcomes that limitation is longduration superpressure balloons (SPBs), which provide observations in a quasiLagrangian frame of reference (e.g., Hertzog et al., 2002; Podglajen et al., 2016a). In the present study, we use SPB measurements as a reference to evaluate the representation of Lagrangian (quasiintrinsicfrequency) spectra. The observations were gathered in the lower stratosphere (16–20 km above sea level) during three SPB campaigns coordinated by the French Space Agency (CNES): Vorcore, PreConcordiasi and Concordiasi. Table 1 summarizes the location and timing of the campaigns, as well as the sampling frequency and status regarding data assimilation. Vorcore and Concordiasi took place in the southern polar vortex during austral spring, while PreConcordiasi flights were launched from an equatorial location (the Seychelles) in boreal spring. Whereas data from the later Concordiasi campaign were broadcast on the Global Telecommunication System (GTS), the earlier Vorcore and PreConcordiasi observations are independent datasets which can be used to evaluate reanalyses (Boccara et al., 2008; Podglajen et al., 2014). For that reason, we will focus in Sect. 3 on Vorcore and PreConcordiasi, while Concordiasi will be used to assess the impact of balloon data assimilation in Sect. 4.1. During the three campaigns, a whole set of measurements was performed, including e.g., ozone and particle measurements; for our purpose, however, only the horizontal position of the balloon from the onboard GPS and pressure and temperature from the thermodynamic sensor (TSEN) will be used.
Apart from occasional changes in the mass or volume of the system (such as dropsonde launching and changes in the incoming radiation flux) and balloon inertia acting at periods shorter than 15 min (Vincent and Hertzog, 2014; Podglajen et al., 2016b), SPBs are essentially passively advected on isopycnic (constantdensity) surfaces once they have ascended to their equilibrium level in the stratosphere. This quasiLagrangian behavior renders measurement interpretation relatively simple. Horizontal winds u and v were derived from the successive positions of the balloon, while pressure is directly measured by the TSEN meteorological sensors. In case of slightly uneven time sampling (varying sampling time step), the raw measurements were interpolated linearly onto a regular time grid. A fast Fourier transform algorithm was then applied to the time series, yielding the signal's Fourier transform ($\widehat{u}\left(\widehat{\mathit{\omega}}\right)$, $\widehat{v}\left(\widehat{\mathit{\omega}}\right)$, $\widehat{p}\left(\widehat{\mathit{\omega}}\right)$). Periodograms were then obtained directly from $\left\widehat{u}\right(\widehat{\mathit{\omega}}){}^{\mathrm{2}}$, $\left\widehat{v}\right(\widehat{\mathit{\omega}}){}^{\mathrm{2}}$ and $\left\widehat{p}\right(\widehat{\mathit{\omega}}){}^{\mathrm{2}}$. In practice, we used a variant of Welch's method and estimate spectra by averaging periodograms obtained from 8 d windows with 4 d overlaps. Note that for consistency with, for example, Fritts and Alexander (2003) our definition of the intrinsic frequency Fourier transform is such that the inverse transform reads as follows:
While this sign convention does not affect the periodograms (estimated from squared modulus quantities and therefore insensitive to the phase), it matters as far as the phase of the signal is concerned, notably for the polarization relations relating the Fourier transforms of the different variables.
Combining the spectra $\left\widehat{u}{}^{\mathrm{2}}\right(\widehat{\mathit{\omega}})$ and $\left\widehat{v}{}^{\mathrm{2}}\right(\widehat{\mathit{\omega}})$ leads to the horizontal kinetic energy spectrum ${E}_{{k}_{\mathrm{h}}}\left(\widehat{\mathit{\omega}}\right)$ per unit mass:
On the other hand, pressure fluctuations p^{′} along the balloon trajectory can be used to estimate the vertical displacement of isopycnic surfaces ${\mathit{\zeta}}_{\mathit{\rho}}^{\prime}$ and isentropic surfaces ${\mathit{\zeta}}_{\mathit{\theta}}^{\prime}$ along airparcel trajectories (see Hertzog et al., 2002; Podglajen et al., 2014, 2017):
where $\stackrel{\mathrm{\u203e}}{\mathit{\rho}}$ is the segmentaveraged density and $\mathit{\alpha}=\left(\frac{g}{{C}_{p}}+\frac{\mathrm{d}\overline{T}}{\mathrm{d}z}\right)/\left(\frac{g}{R}+\frac{\mathrm{d}\overline{T}}{\mathrm{d}z}\right)$ depends on the local temperature lapse rate $\mathrm{d}\overline{T}/\mathrm{d}z$ (g is the gravitational acceleration, C_{p} the thermal capacity of air at constant pressure and R the gas constant for air). The local temperature lapse rate $\mathrm{d}\overline{T}/\mathrm{d}z$ is not directly measured but obtained from the ECMWF ERA5 reanalysis; the sensitivity to the exact value of $\frac{\mathrm{d}\overline{T}}{\mathrm{d}z}$ is small and does not affect the conclusions presented below. Equation (3) relies on a few assumptions (small vertical excursions, small Eulerian pressure perturbation relative to temperature perturbations, adiabatic and hydrostatic flow) which are well met for the low and intermediatefrequency (periods between 1 d and 15 min) motions of interest here (as they are the ones resolvable by the reanalyses). Using Eq. (3), the potential energy spectrum per unit mass E_{p} can be deduced from the estimated pressure spectrum:
As an illustration, Fig. 1 (updated from Podglajen et al., 2016a) shows the intrinsic frequency spectra (power spectral densities or PSDs) of ${E}_{{k}_{\mathrm{h}}}$, E_{p} and of the zonal momentum flux per unit mass $\stackrel{\mathrm{\u203e}}{{u}^{\prime}{w}^{\prime}}$ for the Concordiasi (polar) and PreConcordiasi (tropical) campaigns, which have the highest sampling frequency (Table 1). The important characteristics that are present in the observations and will be searched for in the reanalyses are as follows:

in the polar region, a spectral gap at f between low frequencies and GW frequencies seen in both E_{p} and ${E}_{{k}_{\mathrm{h}}}$ and then a local peak in ${E}_{{k}_{\mathrm{h}}}$ at frequencies just higher than f, almost absent in E_{p}, and in that frequency range ($f<\widehat{\mathit{\omega}}<\mathrm{4}\phantom{\rule{0.125em}{0ex}}f$), ${E}_{\mathrm{p}}<{E}_{{k}_{\mathrm{h}}}$;

in the tropics and in the midfrequency range ($\widehat{\mathit{\omega}}\gg f$) for the polar latitudes, the spectrum follows a ${\widehat{\mathit{\omega}}}^{s}$ power law with s≃2 for E_{p} and ${E}_{{k}_{\mathrm{h}}}$ and s≃1 for $\left{u}^{\prime}{w}^{\prime}\right$, whose scaling appears until $\widehat{\mathit{\omega}}\simeq $ 100 cycles per day (cy d^{−1}).
Above about 100 cy d^{−1}, the balloonobserved E_{p} and $\left{u}^{\prime}{w}^{\prime}\right$ increase and peak near the Brunt–Väisälä frequency N. Although potential physical reasons exists for a spectral peak of momentum and potential energy near N in the atmosphere (Podglajen et al., 2016a), part of the balloonobserved enhancement is likely an artifact caused by the nonisopycnic response of the platform (Podglajen et al., 2016a; Vincent and Hertzog, 2014). Furthermore, due to their expected small horizontal scale and the importance of nonhydrostatic effects (absent from hydrostatic weather models), we can safely assume that buoyancy waves (near N) are absent from hydrostatic, lowresolution weather forecast models. Hence, we will focus our analysis on intrinsic frequencies smaller than 48 cy d^{−1} (periods longer than 30 min). For pressure observations, there is a significant aliasing at 48 cy d^{−1} from higherfrequency motions (in particular near the N peak); this is overcome by computing 15 min Hann windowweighted averages of the highfrequency pressure measurements (recorded every 1 min for Vorcore and 30 s for Concordiasi and PreConcordiasi; see Table 1) and by using this subsampled lowpassfiltered version of the data in the subsequent analysis.
2.2 Atmospheric (re)analyses
For this study, we initially considered four recent reanalysis systems (ERAInterim, ERA5, MERRA2 and JRA55) and the ECMWF operational analysis (corresponding to the model versions used in February and December 2010, i.e., at the times of the PreConcordiasi and Concordiasi balloon campaigns). The model version and spatial resolution of the different products are summarized in Table 2, together with the time resolution at which the outputs (reanalyses and forecasts) are saved and/or made available. For more details regarding the model setups, physical parameterization and data assimilation flow, we refer the reader to the SRIP introduction paper (Fujiwara et al., 2017).
Due to finite storage ability, the (output) time sampling is rather coarse (hourly at best). This may have been a limiting factor for our study. GWs indeed have intrinsic periods ranging from 12 h down to a few minutes (see Fig. 1); so with field outputs every 3 or 6 h, the dominant fraction of GW variance is aliased towards lower frequencies, thus potentially affecting our estimates. However, this limitation is mitigated by the fact that we analyze spectra as a function of intrinsic frequency $\widehat{\mathit{\omega}}$:
which combines the groundbased frequency of the motion ω and its horizontal scale (given by the zonal and meridional wavenumbers k and l). Investigations of different time sampling with ERA5 (for which hourly outputs are available) demonstrated that, while in polar regions the considered intrinsic frequency spectra are strongly sensitive to a change of the sampling time step from 6 to 3 h, using even more frequent (than 3hourly) reanalysis outputs only marginally affects the results. In light of this, the main body of the paper will not present the results obtained with JRA55, available only every 6 h (see Table 2) because of which our spectral analysis suffers from aliasing. For a fair comparison, all (re)analyses (including ERA5) will be used at a time resolution of 3 h except in Appendix A where the impact of the output frequency is investigated.
Gelaro et al. (2017)Kobayashi et al. (2015)Dee et al. (2011)Hersbach et al. (2020)2.3 Comparison methodology
A simplistic approach to compare balloon observations and reanalyses would be to interpolate the model fields along the actual balloon trajectory. However, this might lead to erroneous conclusions regarding the variability present: if, for instance, the balloon were to record a vertical oscillation with a sheared flow u in the reanalysis, this may lead to an oscillation in the interpolated u wind although it might be absent from an actual trajectory computed with the reanalysis wind. To avoid that complication, we directly computed isopycnic (balloonlike) trajectories using the reanalysis fields. In other words, we solve the following system of ordinary differential equations:
in which the air density ρ is a strictly decreasing function of geometric altitude. In practice, System (6) is solved with ln (ρ) as the vertical coordinate, and the 2D trajectories are integrated using a Runge–Kutta method of the order of 4 and a time step of 1 min that is adjusted when needed to satisfy the Courant–Friedrichs–Lewy (CFL) criterion. We note that the dependency of the trajectory on the integration time step and on the details of the numerical scheme is small compared to other sources of errors (Bowman et al., 2013). Gridded reanalysis fields (typically p, T, u and v) are interpolated in the horizontal and time dimensions using cubic splines, leading to vertical profiles. Then the vertical coordinate ln (ρ) is calculated, and finally the wind is interpolated at the density level of the balloon. To examine the Lagrangian variability and estimate spectra, we calculate 8 d trajectories starting at the balloon position every 4 d, thus matching the segments used in the observations. Examples of such trajectories for ERA5 reanalysis are displayed in Fig. 2.
At this point, two remarks should be made. The first is with regards to trajectory accuracy: while it is affected by the sampling frequency (in space and time) of the model and the interpolation method used (Stohl et al., 1995; Bowman et al., 2013), the main source of uncertainty in the lower stratosphere stems from errors in the reanalysis fields (Boccara et al., 2008; Podglajen et al., 2014). In polar regions, the analyses compare well to observations at lower stratospheric altitudes (Boccara et al., 2008), and so there is generally limited difference between observed and simulated balloon trajectories over periods of ∼8 d, as illustrated in Fig. 2 for the case of ERA5. In the tropics, in contrast, analyses may exhibit large deficiencies (Podglajen et al., 2014; Dharmalingam et al., 2019), and the computed trajectories can largely diverge from the observed ones. However, our results are not overly sensitive to the exact number and location of trajectory segments used. Errors in largescale wind are deemed very unlikely to generate a sampling bias and artificially degrade GW variability (which they could do in particular and improbable cases, e.g., by making the balloons drift systematically in quieter areas than those where they actually flew).
Second, we note that physical applications and process studies are interested in air parcel trajectories, which are a priori distinct from the isopycnic trajectories of the balloons. However, as argued above and in Podglajen et al. (2016b), the horizontal positions of isopycnic, isentropic and true air parcel trajectories should remain close (relative to the distance traveled) at synoptic timescales (a few days), and Eq. (4) is expected to be a good approximation in the stratosphere with low diabatic heating. In order to verify this, we compared the 8 d trajectories and spectra for isopycnic, isentropic (constant potential temperature θ) and full 3D trajectories, computed either in p coordinate with $\mathit{\omega}=\mathrm{D}p/\mathrm{D}t$ as vertical velocity (kinematic trajectories) or with θ coordinate and $\dot{\mathit{\theta}}$ vertical velocity (diabatic trajectories). The different horizontal and vertical trajectories for the special case of Vorcore Balloon 2 on 2 November 2005 are displayed in Fig. 2.
Figure 2 demonstrates the validity of our approximations for ERA5: the three trajectory types only slightly diverge from one another in the horizontal. Although as expected the isopycnic trajectory diverges more rapidly from the 3D trajectory than the isentropic one, for an integration time of 8 d the paths remain close. While their vertical positions are clearly distinct, the isentropic, 3D kinematic and isopycnic (scaled onto the isentropes using Eq. 3) trajectories are highly correlated at short timescales. As examined in Appendix B, the estimated Lagrangian energy spectra exhibit only slight differences. Hence, for a more direct comparison with balloon observations, we considered isopycnic trajectories (adjusted using Eq. 3) in the following.
In Fig. 3, the intrinsic frequency spectra (power spectral density) of ${E}_{{k}_{\mathrm{h}}}$, E_{p} and the zonal pseudomomentum flux $(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}})\left{u}^{\prime}{w}^{\prime}\right=(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}})\left\mathrm{\Re}\left(\widehat{u}{\widehat{w}}^{*}\right)\right$ derived from isopycnic trajectories computed using 3hourly wind data from the different (re)analyses are depicted. The trajectories are output every 15 min, so, from the point of view of the sampling, the highest resolvable frequency (Nyquist frequency) in the spectra is 1∕48 cy d^{−1}. However, we should recall that, as mentioned above and explained in Appendix A, the effective resolution in terms of intrinsic frequency actually depends on the limited time and the space resolution of the (re)analyses.
3.1 Horizontal kinetic energy spectra ${E}_{{k}_{\mathrm{h}}}$ $\left(\widehat{\mathit{\omega}}\right)$
Figure 3a and b show the ${E}_{{k}_{\mathrm{h}}}$ spectra (solid lines) obtained from polar (Vorcore; panel a) and equatorial (PreConcordiasi; panel b) trajectories for the balloon observations (black lines) and the reanalyses (colors). In the polar case, lowfrequency features (below the Coriolis frequency f) are wellrepresented in all reanalyses. However, in the gravity wave frequency range (above f), the variance in ${E}_{{k}_{\mathrm{h}}}$ is largely underestimated in all three systems compared to observations. Nevertheless, the reanalysis ${E}_{{k}_{\mathrm{h}}}$ spectra exhibit typical features similar to the observations, including the following:

a local minimum in ${E}_{{k}_{\mathrm{h}}}$ between f∕2 and f, the characteristic spectral gap separating GW ($\widehat{\mathit{\omega}}>f$) from synopticscale motions ($\widehat{\mathit{\omega}}<f$);

a local maximum in ${E}_{{k}_{\mathrm{h}}}$ around 1.2–1.5 f (the nearinertial peak; see Hertzog et al., 2002);

a powerlaw decrease in ${E}_{{k}_{\mathrm{h}}}$ for frequencies larger than 1.5 f.
The differences between reanalyses and observations mainly concern (1) the frequency at which the peak around f in ${E}_{{k}_{\mathrm{h}}}\left(\widehat{\mathit{\omega}}\right)$ occurs and its magnitude and (2) the powerlaw slope of the decreasing ${E}_{{k}_{\mathrm{h}}}$ for $\widehat{\mathit{\omega}}\gg f$; the observations have the shallower slope, followed by ERA5 and ERAInterim, and MERRA2 has the steepest slope.
Despite the quantitative disagreement between the different reanalyses, they also exhibit qualitative similarities with one another and with the observations. In particular, they all show a nearinertial ${E}_{{k}_{\mathrm{h}}}$ peak, which is consistent with observations (Hertzog et al., 2002). In order to further investigate the agreement between the resolved perturbations and gravity wave theory, we examine to what extent the polarization relations for inertiogravity waves are fulfilled. On the f plane, with the convention given by Eq. (1), the polarization relation for the horizontal wind reads as follows (Fritts and Alexander, 2003):
where ${\widehat{u}}_{\parallel}$ is the amplitude of the horizontal wind along the wave vector and ${\widehat{u}}_{\u27c2}$ the amplitude perpendicular to the wave vector. Keeping in mind the convention expressed by Eq. (1), Eq. (7) indicates that lowfrequency waves (near f) induce anticyclonic flow rotation, whereas highfrequency waves have their horizontal wind perturbation aligned with the wave vector and no preferred rotation direction. To make use of this property and further demonstrate that the nearinertial peak in the analyses is due to inertial oscillations, we turn to rotary spectral analysis. The rotary Fourier transform of the horizontal wind $\widehat{U}$ is defined by
As a consequence of Eq. (7), the rotary ratio $R\left(\widehat{\mathit{\omega}}\right)$, which is the ratio of the PSD of counterclockwise horizontal motions, $\left\widehat{U}\right(\widehat{\mathit{\omega}}){}^{\mathrm{2}}$ with our Fourier transform convention in Eq. (1), and clockwise ones ($\left\widehat{U}\right(\widehat{\mathit{\omega}}){}^{\mathrm{2}}$), verifies the relation:
Note that the rotary ratio depends on the sign of f, so the rotation direction is reversed in the Southern Hemisphere with respect to the Northern Hemisphere, with anticyclonic motions being always favored.
Equation (9) has been exploited in previous studies (e.g., Hertzog et al., 2002; Conway et al., 2019) to evaluate the consistency of the observed horizontal wind spectrum with linear gravity wave theory. In particular, the dominance of anticyclonic (here in the Southern Hemisphere) motions at low frequencies is consistent with Eq. (9) and indicates the importance of linear gravity waves as opposed to, for example, stratified turbulence for which no systematic phase relation between the horizontal wind components is expected. In Fig. 3, $R\left(\widehat{\mathit{\omega}}\right)$ is displayed together with the theoretical value ratio given by Eq. (9). The dominance of anticyclonic motions is a further argument in favor of inertioGWs being responsible for the f peak in both reanalyses and observations. However, we note that it is less pronounced in the reanalyses than in the observations and in ERA5 than ERAInterim and MERRA2. Possible reasons for this will be examined in the next section.
In the tropics, the statistical agreement between reanalyses and observations regarding the representation of highfrequency variability is surprisingly better than in the polar latitudes. This contrasts with the representation of largescale wind fields which is better in polar regions (Boccara et al., 2008) than in equatorial ones (Podglajen et al., 2014). The agreement obtained for lowfrequency waves is quantitative as well as qualitative and reaches up to a cutoff frequency which is about 4 cy d^{−1} for the ERAInterim, ERA5 and the ECMWF operational analyses and 2 cy d^{−1} for MERRA2. Below this cutoff frequency, the observed and modeled age spectra are in quantitative and qualitative agreement. At the cutoff frequency, the reanalysis spectra show a kink, and above it they drop with much larger slopes than the observed ones, which keep decaying with a constant slope.
3.2 Potential energy spectra
Figure 3c and d show the spectra of the potential energy per unit mass E_{p} for the polar (panel c) and tropical (panel d) flights. To a large extent, the situation is similar to the one described above for ${E}_{{k}_{\mathrm{h}}}$: in the polar case, there is a qualitative agreement in the structure of the power spectra (notably regarding the spectral gap around f) but a quantitative disagreement and a steeper powerlaw slope than in the observations.
Besides the direct value of the potential energy, it is again interesting to investigate whether the fluctuations in reanalyses verify the polarization relation expected for gravity waves. In the absence of constructive interference, the ratio of potential to horizontal kinetic energy for frequencies $\widehat{\mathit{\omega}}\ll N$ should obey (Podglajen et al., 2016b)
Figure 3 exhibits $\frac{{E}_{{k}_{\mathrm{h}}}\left(\widehat{\mathit{\omega}}\right)}{{E}_{\mathrm{p}}\left(\widehat{\mathit{\omega}}\right)}$ for the three reanalyses and the observations (dashed lines), together with its theoretical value given in Eq. (10). The observations again closely follow theoretical expectations from ∼1.2 f, with a dominance of kinetic energy, to higher frequencies (the socalled midfrequency range) at which there is an equipartition between E_{p} and ${E}_{{k}_{\mathrm{h}}}$. The spectra of the fluctuations along reanalysis trajectories show the same evolution with frequency as the observations and the ${E}_{{k}_{\mathrm{h}}}/{E}_{\mathrm{p}}$ ratio decreasing with increasing frequency. However, reanalyses also suffer from a significant overestimation of the kinetic energy compared to the potential energy. This is the case by a factor of 1.5 to 3 for ERA5 and 5 for ERAInterim and MERRA2. Together with the rotary ratio analysis, this suggests that a significant fraction of the variability does not obey the observed (and expected) polarization relations for gravity waves. One might speculate that this inconsistency stems from numerical dissipation.
In the tropics, as for the ${E}_{{k}_{\mathrm{h}}}$ spectra, the E_{p} spectra from reanalyses and observations are in better agreement than over the pole. In particular, the ${E}_{{k}_{\mathrm{h}}}/{E}_{\mathrm{p}}$ ratio is close to 1 for the reanalyses, as expected, whereas there is a quantitative agreement between observed and analyzed E_{p} up to a threshold frequency similar to that encountered with ${E}_{{k}_{\mathrm{h}}}$. Besides the reanalyses, for the tropical case the spectra from the ECMWF operational analysis trajectories are also displayed; they show similar quantitative results as for the ERA5 reanalysis.
3.3 Eliassen–Palm flux and vertical wind spectra
Figure 3e and f display the zonal pseudomomentum $\left(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}}\right)\left{u}^{\prime}{w}^{\prime}\right$ flux spectra of resolved waves in the reanalyses. This quantity characterizes the forcing of the largescale flow by the waves, so these panels enlighten us regarding the missing gravity wave drag from resolved waves compared to observations.
Over the pole, the reanalyses underestimate the variability in the whole GW range, with ERA5 being the closest to the balloon observations. Although the number of reanalysis products is limited, there appears to be a correlation between the fraction of variability resolved at low frequency and the vertical resolution of the product (given in Table 2). With the highest vertical resolution, ERA5 resolved the largest fraction of the variability, followed by ERAInterim and MERRA2 with coarser resolutions.
In the equatorial region, on the other hand, similar to the ${E}_{{k}_{\mathrm{h}}}$ and E_{p} spectra, the momentum flux by lowfrequency equatorial GWs is comparable to observations in all products up to the threshold frequency where it drops with a slope twice as large as the actual slope.
We do not show the vertical wind spectrum, but its shape can be readily deduced from the potential energy spectrum as ${E}_{{k}_{\mathrm{v}}}=\frac{\mathrm{1}}{\mathrm{2}}{\widehat{w}}^{\mathrm{2}}\propto {\left(\frac{\widehat{\mathit{\omega}}}{N}\right)}^{\mathrm{2}}{E}_{\mathrm{p}}$, so the unequal performances of the reanalyses can be deduced from Fig. 3c and d. It should be noted that the different quantities considered have various powerlaw slopes in the gravity wave range from ${\widehat{\mathit{\omega}}}^{\mathrm{2}}$ for the horizontal kinetic energy ${E}_{{k}_{\mathrm{h}}}$ and potential energy E_{p} power spectra, to ${\widehat{\mathit{\omega}}}^{\mathrm{1}}$ for the Eliassen–Palm (EP) flux spectrum $\left(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}}\right)\left{u}^{\prime}{w}^{\prime}\right$, and to ${\widehat{\mathit{\omega}}}^{\sim \mathrm{0}}$ for the vertical kinetic energy ${E}_{{k}_{\mathrm{v}}}$. As a consequence, the variability of different fields emphasizes different parts of the spectrum; while u, v (${E}_{{k}_{\mathrm{h}}}$) and T (E_{p}) are more connected with the lowfrequency part of the gravity wave spectrum, $\left(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}}\right)\left{u}^{\prime}{w}^{\prime}\right$ corresponds to the intermediate frequencies, and the vertical wind component w (related to ${E}_{{k}_{\mathrm{v}}}$) corresponds to the highfrequency waves.
3.4 Intermittency and distribution of the fluctuations
The PSDs examined above provide information on the autocorrelation and the repartition of the fluctuations as a function of intrinsic frequency. However, as pointed out by several authors (e.g., Hertzog et al., 2012; Podglajen et al., 2016b), the probability distributions are also relevant to determine whether the variance is due to ubiquitous, constant variability generated by the superposition of random waves or if rare, large excursions created by specific wave events can occur. The latter behavior is expected for GWs which are known to be intermittent.
3.4.1 Intermittency of the fluctuations
Similar to what was used for the Lagrangian spectra, we compare the distribution of the fluctuations in the reanalyses to the observations based on computed isopycnic trajectories corrected for the nonisopycnic behavior (by using the coefficient α defined in Eq. 3). In order to focus on the GW range, we filter the outputs of the trajectories to keep only the signal corresponding to intrinsic frequencies shorter than 48 cy d^{−1} and larger than $f\left(\mathrm{30}{}^{\circ}\right)/\left(\mathrm{2}\mathit{\pi}\right)\simeq \mathrm{1}$ cy d^{−1} (tropics) or $f\left(\mathrm{70}{}^{\circ}\right)/\left(\mathrm{2}\mathit{\pi}\right)\simeq \mathrm{2}$ cy d^{−1} (pole). The lower bound is considered the lower bound of the GW frequency range in the regions studied, while the upper bound corresponds to the Nyquist frequency of Vorcore balloon position information (every 15 min). It is below the frequency at which the balloons start to depart significantly from the isopycnic behavior, so the comparison can only be marginally affected by the nonisopycnic balloon response. In that frequency range $[\frac{f}{\mathrm{2}\mathit{\pi}};\sim \frac{N}{\mathrm{4}\mathit{\pi}}]$, we consider temperature, zonal momentum flux and vertical velocity to characterize statistics and intermittency of low, medium and highfrequency waves, respectively.
Figure 4 shows the probability distribution (PDF) of the three quantities (temperature, zonal momentum flux and vertical wind) for all the (re)analyses considered. From top to bottom, the considered fluctuations are primarily induced by waves of higher frequency; ${T}_{l}^{\prime}{}^{\mathrm{2}}\propto {E}_{\mathrm{p}}(\widehat{\mathit{\omega}})$ emphasizes the low frequencies, while $\left{u}^{\prime}{w}^{\prime}\right\propto \widehat{\mathit{\omega}}{E}_{\mathrm{p}}\left(\widehat{\mathit{\omega}}\right)$ and ${w}^{\prime}{}^{\mathrm{2}}\propto {\widehat{\mathit{\omega}}}^{\mathrm{2}}{E}_{\mathrm{p}}(\widehat{\mathit{\omega}})$ are related to increasing frequencies. Since highfrequency waves are more poorly represented than lowfrequency ones (Fig. 3), it comes as no surprise that the PDF of analyzed and observed fluctuations is in increasing disagreement from temperature to momentum flux and from momentum flux to vertical velocity in both tropical and southern polar regions.
The width of the temperature PDFs for the polar and tropical regions is redundant with the level of the variance in the E_{p} spectrum already shown in Fig. 3. The additional information provided in Fig. 4 concerns the shape of the PDF. In the tropics, GW temperature fluctuations are ubiquitous and characterized by a Gaussian PDF. This behavior is fairly well reproduced in all reanalyses examined here. In contrast, at the pole, the balloon PDFs are no longer Gaussian. In particular, they have longer tails than Gaussian PDFs, a behavior which is reproduced by reanalyses.
Regarding the momentum flux, the tropical and polar PDFs both exhibit lognormal PDFs, as well as the reanalyses. We note that the tropical regions during PreConcordiasi show more activity than the polar regions during Vorcore, but this is largely due to a wider GW range in the tropics, where f→0. Finally, the vertical wind PDFs show Laplace distributions in both cases, and the order of the reanalyses is the same as for the rest. In summary, while the fraction of resolved variance varies from one analysis to the other, it appears that the basic intermittency properties of the GW field and the shape of the PDF of the fluctuations are consistent between observations and the different reanalyses.
It should also be noted that the conclusion regarding the resolved fraction of variability for the different reanalyses is transferable from one campaign to another. In other terms, for all quantities and regions examined here, we find that the more realistic reanalysis is ERA5, followed by ERAInterim and MERRA2. The 2010 ECMWF operational model performs better than ERA5 likely because of its superior horizontal resolution.
3.4.2 Geographic distribution of the fluctuations
We have considered above the intermittency of GWs sampled by the balloons along their trajectories. This intermittency stems from both the space and time variability of the GW field in which the balloons drift. A significant part of it comes from the geographic variability of wave sources, as illustrated in Fig. 5 through global maps of the vertical wind standard deviation at 100 hPa (related to highfrequency wave activity) for the month of March 2010. The geographic structure is similar between the different reanalyses, which is consistent with the fact that intermittency is similar in the different products (see the shapes in Fig. 4). GW amplitudes largely differ from one product to the other, but the geographic variability is to a large extent common to all reanalyses. In the season examined (boreal spring 2010), mountain ranges such as the Rockies, the Andes, the Himalayas and the Antarctic Peninsula stand out as regions of increased activity. Convective regions such as the Intertropical Convergence Zone are also characterized by a larger activity in all products. However, although the general geographic pattern of σ_{w} does match the different products, there are differences in the details. In particular, the features are sharper and smaller scale (especially around orography) in the highresolution products (the operational analysis and ERA5) than in the lowerresolution reanalyses (ERAInterim, MERRA2 and JRA55). Although this property might be expected, it implies that intermittency is increased in those products.
Since the different resolutions in the reanalyses imply different horizontal wave numbers which may undergo different filtering and vertical propagation, we present a vertical profile of σ_{w} in the tropics for the different products in Fig. 6a. There is a drastic reduction in σ_{w} from the troposphere to the stratosphere related to the increased stratospheric stability hampering vertical motions there. Figure 6b displays the same σ_{w} equatorial profile but normalized by the standard deviation from the ECMWF analysis. Although they do not resolve the same wave population, the different analyses essentially show a similar vertical structure in GW variance except ERA5 for which σ_{w} is relatively reduced around the tropopause. This is likely due to the increased stability in that region in the high vertical resolution ERA5, which can better resolve strong vertical gradients near the tropopause.
4.1 Impact of balloon data assimilation
Balloon observations from the PreConcordiasi and Vorcore campaigns were not assimilated in the reanalyses, so they provide an independent evaluation of the resolved GW variability. However, data collected during the Concordiasi campaign, which took place in austral winter 2010, were broadcast on the GTS and assimilated in most analyses, including ECMWF operational data, ERA5, ERAInterim and MERRA2 (Hoffmann et al., 2017). Hence, comparing the period of the (assimilated) Concordiasi campaign with that of the (unassimilated) Vorcore campaign provides an opportunity to characterize the extent to which balloon data assimilation in the atmospheric models impacts the representation of the wave field (either through spurious wave generation or realistic waves introduced in the initial state that propagate in the forecast). Since Concordiasi only flew in austral spring 2010, GW activity in the 2010 southern lower stratospheric polar vortex is not typical but largely influenced by data assimilation. However, with the development of the Loon dataset of superpressure balloons (Conway et al., 2019) and studies regarding its assimilation in numerical weather prediction (NWP) models (Coy et al., 2019), it is interesting to document the impact of such a SPB dataset.
The ${E}_{{k}_{\mathrm{h}}}$ spectra for the two campaigns are shown in Fig. 7. There are slight differences between the two observed spectra (black lines), notably with a more pronounced f peak during Vorcore than Concordiasi. Those are likely due to the different latitudinal sampling during the two campaigns, with more variable f for Vorcore trajectories than Concordiasi's (see the shaded gray area in Fig. 7). However, the main differences between the two campaigns do not lie in the observed spectra but rather in the analyzed ones. Whereas the reanalysis ${E}_{{k}_{\mathrm{h}}}$ spectra during Vorcore largely underestimate ${E}_{{k}_{\mathrm{h}}}$ and in particular the magnitude of the spectral peak near f, there is a clearer spectral peak in the reanalyses with assimilated balloon data. This better performance of the reanalyses during Concordiasi shows that, besides adjusting the state of the flow, assimilation of the balloon dataset statistically reinforces some dynamical ingredients in the simulation (here inertial waves).
4.2 Impact of underlying model version and resolution on GW representation
4.2.1 Time sampling
Proper representation of the GW variability requires sufficient time resolution particularly to avoid aliasing of the spectral gap and the f peak. As mentioned above, this is not as much of a critical issue for the model time step (which is small enough to describe the temporal evolution of the spatially resolved GW) as for its time sampling. The time sampling of the reanalysis should indeed be high enough to isolate unequivocally the highest groundbased frequency present in the simulation. Furthermore, this parameter is one of the most sensitive in trajectory calculations (e.g., Pisso et al., 2010; Bowman et al., 2013) at least with the spatial resolution and integration time step we used. This is especially the case since the winds are instantaneous fields rather than time averages, which implies aliasing when subsampling (e.g., Stohl et al., 1995; Bowman et al., 2013). In order to test the impact of this key parameter on our Lagrangian spectrum estimation, we performed trajectory calculations with ERA5 varying the time sampling of the reanalysis outputs between 6, 3 and 1 h. We use the instantaneous ERA5 wind values and do not apply any time averaging which would reduce aliasing but which cannot be applied to products with coarser temporal sampling (Hoffmann et al., 2019).
The results of this experiment are presented in Appendix A and can be summarized as follows. First, improving the time sampling of the reanalysis from every 6 h to every 3 h induces large changes in the Lagrangian GW spectrum over the pole and enables inertial oscillations at periods of ∼12 h to be resolved. In contrast, a further improvement of the time resolution from every 3 h to every 1 h has little noticeable effect, suggesting that a large fraction of the resolved high intrinsic frequency GW activity is actually caused by the background wind moving the air parcels across “quasistationary” wave features with groundbased periods longer than about 6 h.
4.2.2 Impact of vertical and horizontal resolution
The nominal resolution of the model, given in Table 2, largely controls the magnitude of GW fluctuations in reanalyses. This is shown in Fig. 8, which displays the fraction of resolved variance in the 2–48 cy d^{−1} intrinsic frequency range in the reanalysis products compared to balloon observations for the horizontal kinetic energy ${E}_{{k}_{\mathrm{h}}}$, the potential energy E_{p}, the zonal pseudomomentum flux $\left(\mathrm{1}{f}^{\mathrm{2}}/{\widehat{\mathit{\omega}}}^{\mathrm{2}}\right)\left{u}^{\prime}{w}^{\prime}\right$ and the vertical kinetic energy ${E}_{{k}_{\mathrm{v}}}$. Considering the (re)analyses from ECMWF (i.e., ERAInterim, ERA5 and op. ECMWF), the sorting of the resolved variance reflects the horizontal resolution of the products for all considered variables. MERRA2 stands out in this consideration; although its nominal resolution (∼60 km; see Table 2) is better than that of ERAInterim (∼79 km), it systematically resolves a smaller fraction of GW variance. However, this reanalysis has a nonspectral grid, which likely implies some additional diffusion to ensure numerical stability, thus reducing the effective resolution of the model (Holt et al., 2016).
Although the dependency on the reanalysis horizontal resolution is present for all variables considered, it is more pronounced at the lowest intrinsic frequency (i.e., for the variables on the left in Fig. 8). Indeed, the dependency on horizontal resolution appears to be stronger for ${E}_{{k}_{\mathrm{h}}}$ than for variables with variance primarily contained at large $\widehat{\mathit{\omega}}$ ($\stackrel{\mathrm{\u203e}}{{u}^{\prime}{w}^{\prime}}$ and ${E}_{{k}_{\mathrm{v}}}$). This is expected when acknowledging that $\widehat{\mathit{\omega}}$ is loosely related to the horizontal wavenumber; while ${E}_{{k}_{\mathrm{h}}}$ and E_{p} are already partly resolved in lowresolution products and do not depend as much on resolution, ${E}_{{k}_{\mathrm{v}}}$ increases strongly when including the additional small scales brought by the high resolution. Along the same line of argumentation, we note that the fraction of resolved variance is larger for the variables depending on lowfrequency waves (${E}_{{k}_{\mathrm{h}}}$ and E_{p}) compared to those whose variance is focused at high frequency. In particular, the fraction of resolved E_{p} is closer to the fraction of resolved ${E}_{{k}_{\mathrm{h}}}$ in the tropics than over the pole, where there is a relative increase in ${E}_{{k}_{\mathrm{h}}}$ variance at low frequency which is better resolved than the variability at a higher frequency. The difference seen between the tropical and polar region in terms of resolved ${E}_{{k}_{\mathrm{h}}}$ and E_{p} variance is almost absent for ${E}_{{k}_{\mathrm{v}}}$.
Since GWs occur at the mesoscale where the gray zone of global models currently lies, the strong sensitivity to horizontal resolution is expected; models with higher resolution are able to resolve a larger fraction of the spectrum due to their higher cutoff wavenumber. For a fair comparison of the datasets with different resolutions, we also considered a version of ERA5 truncated in spectral space to match the resolution of ERAInterim (truncation T255). The GWs present are those that can be represented on the ERAInterim grid; hence, differences arise from the different vertical resolutions, propagation properties and sources of largescale horizontal waves. Note that, since it still plays a role in the generation and propagation of the resolved waves, this does not remove the impact of horizontal resolution entirely but rather provides a lower bound estimate of this impact. The corresponding map of vertical wind standard deviation ${\mathit{\sigma}}_{\mathrm{w}}=\sqrt{\mathrm{2}\phantom{\rule{0.125em}{0ex}}{E}_{{k}_{\mathrm{v}}}}$ is shown in Fig. 5d. Two interesting results emerge from this exercise. First, the confinement of regions of highGW activity is similar between closetonative resolution ERA5 and its spectrally truncated version but much less pronounced in ERAInterim where regions of large GW activity are more spread around. Second, although they effectively have the same resolution, the truncated ERA5 has enhanced variance compared to ERAInterim. Besides horizontal resolution, vertical resolution and vertical mixing formulation have also been shown to play a dominant role in controlling the wind energetics in atmospheric models (Skamarock et al., 2019). In particular, in their global simulations, Skamarock et al. (2019) found that both the horizontal kinetic energy spectrum and the wind shear (quantified by the gradient Richardson number) largely increased in amplitude with vertical resolution, with the first signs of convergence observed with vertical mesh spacing near 100 m. Nearinertial waves with reduced vertical scales are believed to contribute greatly to the enhanced wind variance at high vertical resolutions (Waite and Snyder, 2009; Skamarock et al., 2019). The dispersion relation,
indeed shows that when $\widehat{\mathit{\omega}}\to f$, the vertical wavelength is reduced, and these waves are potentially more sensitive to the vertical resolution. In accordance with those previous works, our analysis suggests a specific difficulty in simulating small vertical wavelength waves near the inertial frequency. The ${E}_{{k}_{\mathrm{h}}}$ amplitude of GWs of a frequency between 2 and 4 cy d^{−1} is only slightly underestimated in the tropics (see Fig. 3) where they have larger vertical scales (f→0 in Eq. 11), whereas it is much reduced in the polar latitude (Fig. 3) where those waves are nearinertial waves levels and hence close to critical levels. However, the different reanalysis products from the ECMWF do not exhibit large changes in the description of the f peak although the vertical resolution of the model in the lower stratosphere varies by more than a factor of 4 between ERAInterim and ERA5 (Table 2). This is the case whether or not balloon data are assimilated (Fig. 7). When considering Fig. 8, there is no clear (i.e., monotonic) relationship between vertical resolution and the fraction of resolved variance, with, for instance, the ECMWF operational analysis performing better in the tropics than the higher vertical resolution ERA5. It is likely that the limited vertical resolution of the reanalyses in the lower stratosphere is too low or the vertical mixing too large for the sensitivity to vertical resolution to clearly emerge.
4.3 Implications for gravity wave studies and Lagrangian modeling based on reanalyses
A simple recommendation can be drawn from the analysis presented above. As pointed out by a number of studies, gravity waves are now partly represented in modern reanalyses, in particular if one considers the temperature and horizontal wind fluctuations which are tied to the lowfrequency part of the gravity wave intrinsic frequency spectrum. In the context of lowfrequency largescale waves, the use of reanalysis to quantify atmospheric gravity wave variability is justified in particular for case studies in the tropics. However, the fraction of the variability resolved by the reanalyses gets smaller with the scale of the waves that generate that variability, so the momentum flux (intermediatescale waves) and a fortiori the vertical wind variability (tied to smallscale waves) are still underestimated.
This paper examined the representation of gravitywaveinduced fluctuations in the lower stratosphere of modern reanalyses through the computation of Lagrangian trajectories and their comparison with quasiLagrangian balloon observations. Consistent with previous studies, we find that reanalysis systems resolve a significant part of the gravity wave spectrum at low intrinsic frequencies and represent part of the associated variability in temperature and horizontal wind. In particular, specific characteristics of the observed spectrum can be well reproduced, such as the spectral gap between gravity waves and the largescale flow or the increased horizontal wind variance near the Coriolis frequency. However, the resolved fraction of the variability decreases with increasing intrinsic frequency so that the vertical wind variability, tied to high intrinsic frequency, is largely underestimated in modern systems. The momentum flux, peaking at intermediate scales, is also underestimated though less than vertical velocity. Generally, the analyses examined fulfill the expectation that the higher the nominal horizontal resolution, the larger the fraction of resolved variability. Contrasting the tropics with the southern high latitudes further suggests that small vertical wavelength waves (approaching the inertial frequency) at high latitudes are not well represented. Despite the capability of the underlying models, we also found that the waves in the analyses are not purely self generated by the models but may rely for a significant part on data assimilation. In particular we found that the assimilation of longduration balloon observations over Antarctica in austral spring enhanced GW activity in the reanalyses. In any case, the fact that GWs are present in NWP models suggests that observations with GW signature (radiosondes, GPS radiooccultations) could be successfully assimilated.
In the future, improved vertical and horizontal resolutions (and maybe specific data assimilation strategies) will allow NWP models to resolve an increased fraction of GWinduced variability, which will make them even more valuable than currently for GWrelated studies. The increased presence of this highfrequency horizontal wind variability requires the use of higher time sampling (3 hourly) for accurate trajectory calculations. However, we expect that vertical wind variability, which is tied to the smallest spatial scales (∼1–10 km) and for which current modeling approximations (such as the hydrostatic one) are not valid, will still remain a challenge. This emphasizes the longterm need for gravity wave parameterizations in NWPbased process studies, especially for processes depending on vertical velocity.
As indicated by Eq. (5), a range of horizontal wavenumbers and ground relative frequencies corresponds to a given intrinsic frequency. It is unclear what the contribution of low and high horizontal phase speed is and whether part of the missing activity results from undersampling the time space. Indeed, the horizontal kinetic energy, for example, at a given intrinsic frequency ${E}_{{k}_{\mathrm{h}}}\left(\widehat{\mathit{\omega}}\right)$ actually depends on different groundbased frequencies:
where δ is the 1D Dirac distribution and E is the 3D power spectral density (note that only 3D integration is necessary since for internal gravity waves the dispersion relation directly imposes the vertical wavenumber given the three other wave parameters and the background). It is apparent that a given $\widehat{\mathit{\omega}}$ corresponds to a whole surface of ω and wavenumbers. High intrinsic frequency waves may thus have a high ground relative frequency or be shortscale waves in a large background wind. The two cases are affected differently by our analysis method. Since we retrieve the reanalysis at a horizontal resolution close to the actual grid spacing, we fully resolve analysis waves which have high intrinsic frequencies due to their small horizontal scales (large $k\overline{u}$). However, the reanalyses are sampled at a time resolution of typically 1 to 6 h, i.e., lower than the model time step (less than 15 min). There can be aliasing in the $\widehat{\mathit{\omega}}$ for highfrequency waves with large scales but high ω. In order to evaluate this effect, we take advantage of the 1hourly ERA5 to compare trajectories using 1, 3 or 6hourly reanalysis fields.
The results are displayed in Fig. A1 for Vorcore flights. Large differences occur for both the ${E}_{{k}_{\mathrm{h}}}$ spectrum and the rotary ratio between 6 and 3 h; the inertial peak at $\widehat{\mathit{\omega}}\simeq \mathrm{12}$ h is completely aliased with the resolution of 6 h (for which it is close to the Nyquist frequency) but resolved at 3 h. This shows that the waves dominating the energy spectrum at frequencies near f have predominantly large scales and $\widehat{\mathit{\omega}}\simeq \mathit{\omega}$. In contrast, further improving does not affect the ${E}_{{k}_{\mathrm{h}}}$ spectrum much but leads to a more realistic rotary ratio at higher frequencies. This exercise demonstrates that the 3hourly resolution provided by most reanalyses (except JRA55) is sufficient for our purpose, i.e., evaluating the GW intrinsic frequency spectrum.
To test the validity of the assumptions made in Sect. 2.1 to use balloon isopycnic trajectories to infer Lagrangian spectra, we computed isopycnic, isentropic, and 3D kinematic and 3D diabatic trajectories in the reanalyses. The resulting ${E}_{{k}_{\mathrm{h}}}$ and E_{p} spectra are presented in Fig. B1 for ERA5 with two time resolutions: 1 and 3 hourly. For the ${E}_{{k}_{\mathrm{h}}}$ spectrum, the four calculations converge at 1 and 3 hourly. For the E_{p} spectrum, only the kinematic case with 3hourly winds disagrees with the others. The discrepancy is resolved at 1hourly resolution, suggesting aliasing of highfrequency vertical motions at 3hourly resolution; this effect is filtered out in the isopycnic, isentropic and diabatic cases, for which vertical wind is implicitly accounted through the fully resolved displacement of the θ surfaces.
The balloon data used in this study are archived by the French Atmospheric Data Center (Aeris) and can be located through their catalog at https://www.aerisdata.fr/catalogue/ (last access: 5 August 2020). They are currently available at ftp://aerisgd:aeris@ftp.ipsl.fr/BALLOON/RUMBA_L2 (Hertzog and Vial, 2020) (Vorcore) and ftp://aerisgd:aeris@ftp.ipsl.fr/BALLOON/TSEN_L2 (Hertzog, 2020) (Concordiasi). The (re)analysis datasets can be retrieved from the respective forecast center.
AP designed the study and carried out the analysis with important contributions from AH and RP. AH and AP developed the balloon trajectory code used in the study. BL guided the use of reanalysis data. AP wrote the paper with contributions from all coauthors.
The authors declare that they have no conflict of interest.
This article is part of the special issue “The SPARC Reanalysis Intercomparison Project (SRIP) (ACP/ESSD interjournal SI)”. It is not associated with a conference.
Aurélien Podglajen thanks Manfred Ern for his insightful comments on an earlier version of the paper. The authors are grateful to Corwin Wright and an anonymous referee for their thorough review which helped clarify the paper. ERA5 results were generated using Copernicus Climate Change Service information 2019.
The article processing charges for this openaccess publication were covered by a research center of the Helmholtz Association.
This paper was edited by Gabriele Stiller and reviewed by Corwin Wright and one anonymous referee.
Andrews, D., Holton, J., and Leovy, C.: Middle Atmosphere Dynamics, International geophysics series, Academic Press, available at: https://books.google.fr/books?id=N1oNurYZefAC (last access: 5 August 1010), 1987. a
Bacmeister, J. T., Eckermann, S. D., Tsias, A., Carslaw, K. S., and Peter, T.: Mesoscale Temperature Fluctuations Induced by a Spectrum of Gravity Waves: A Comparison of Parameterizations and Their Impact on Stratospheric Microphysics, J. Atmos. Sci., 56, 1913–1924, https://doi.org/10.1175/15200469(1999)056<1913:MTFIBA>2.0.CO;2, 1999. a
Baldwin, M. P., Gray, L. J., Dunkerton, T. J., Hamilton, K., Haynes, P. H., Randel, W. J., Holton, J. R., Alexander, M. J., Hirota, I., Horinouchi, T., Jones, D. B. A., Kinnersley, J. S., Marquardt, C., Sato, K., and Takahashi, M.: The quasibiennial oscillation, Rev. Geophys., 39, 179–229, https://doi.org/10.1029/1999RG000073, 2001. a
Boccara, G., Hertzog, A., Basdevant, C., and Vial, F.: Accuracy of NCEP/NCAR reanalyses and ECMWF analyses in the lower stratosphere over Antarctica in 2005, J. Geophys. Res., 113, D20115, https://doi.org/10.1029/2008JD010116, 2008. a, b, c, d
Bowman, K. P., Lin, J. C., Stohl, A., Draxler, R., Konopka, P., Andrews, A., and Brunner, D.: Input Data Requirements for Lagrangian Trajectory Models, B. Am. Meteorol. Soc., 94, 1051–1058, https://doi.org/10.1175/BAMSD1200076.1, 2013. a, b, c, d
Conway, J. P., Bodeker, G. E., Waugh, D. W., Murphy, D. J., Cameron, C., and Lewis, J.: Using Project Loon Superpressure Balloon Observations to Investigate the Inertial Peak in the Intrinsic Wind Spectrum in the Midlatitude Stratosphere, J. Geophys. Res.Atmo., 124, 8594–8604, https://doi.org/10.1029/2018JD030195, 2019. a, b
Coy, L., Schoeberl, M. R., Pawson, S., Candido, S., and Carver, R. W.: Global Assimilation of Loon Stratospheric Balloon Observations, J. Geophys. Res.Atmos., 124, 3005–3019, https://doi.org/10.1029/2018JD029673, 2019. a
Dee, D. P., Uppala, S. M., Simmons, A. J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M. A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., MongeSanz, B. M., Morcrette, J.J., Park, B.K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.N., and Vitart, F.: The ERAInterim reanalysis: configuration and performance of the data assimilation system, Q. J. Roy. Meteor. Soc., 137, 553–597, https://doi.org/10.1002/qj.828, 2011. a
Dharmalingam, S., Plougonven, R., Hertzog, A., Podglajen, A., Rennie, M., Isaksen, L., and Kébir, S.: Accuracy of Balloon Trajectory Forecasts in the Lower Stratosphere, Atmosphere, 10, https://doi.org/10.3390/atmos10020102, 2019. a
Dinh, T., Podglajen, A., Hertzog, A., Legras, B., and Plougonven, R.: Effect of gravity wave temperature fluctuations on homogeneous ice nucleation in the tropical tropopause layer, Atmos. Chem. Phys., 16, 35–46, https://doi.org/10.5194/acp16352016, 2016. a
Fritts, D. C. and Alexander, M. J.: Gravity wave dynamics and effects in the middle atmosphere, Rev. Geophys., 41, 1003, https://doi.org/10.1029/2001RG000106, 2003. a, b
Fujiwara, M., Wright, J. S., Manney, G. L., Gray, L. J., Anstey, J., Birner, T., Davis, S., Gerber, E. P., Harvey, V. L., Hegglin, M. I., Homeyer, C. R., Knox, J. A., Krüger, K., Lambert, A., Long, C. S., Martineau, P., Molod, A., MongeSanz, B. M., Santee, M. L., Tegtmeier, S., Chabrillat, S., Tan, D. G. H., Jackson, D. R., Polavarapu, S., Compo, G. P., Dragani, R., Ebisuzaki, W., Harada, Y., Kobayashi, C., McCarty, W., Onogi, K., Pawson, S., Simmons, A., Wargan, K., Whitaker, J. S., and Zou, C.Z.: Introduction to the SPARC Reanalysis Intercomparison Project (SRIP) and overview of the reanalysis systems, Atmos. Chem. Phys., 17, 1417–1452, https://doi.org/10.5194/acp1714172017, 2017. a, b
Gary, B. L.: Mesoscale temperature fluctuations in the stratosphere, Atmos. Chem. Phys., 6, 4577–4589, https://doi.org/10.5194/acp645772006, 2006. a
Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.: The ModernEra Retrospective Analysis for Research and Applications, Version 2 (MERRA2), J. Climate, 30, 5419–5454, https://doi.org/10.1175/JCLID160758.1, 2017. a
Geller, M. A., Alexander, M. J., Love, P. T., Bacmeister, J., Ern, M., Hertzog, A., Manzini, E., Preusse, P., Sato, K., Scaife, A. A., and Zhou, T.: A Comparison between Gravity Wave Momentum Fluxes in Observations and Climate Models, J. Climate, 26, 6383–6405, https://doi.org/10.1175/JCLID1200545.1, 2013. a
Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., MuñozSabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E. Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999–2049, https://doi.org/10.1002/qj.3803, 2020. a
Hertzog, A.: In situ measurements of atmospheric parameters using TSEN instrument, available atL ftp://aerisgd:aeris@ftp.ipsl.fr/BALLOON/TSEN_L2, last access: 5 August 2020. a
Hertzog, A. and Vial, A.: In situ measurements of atmospheric parameters using Nacelle RUMBA, available at: ftp://aerisgd:aeris@ftp.ipsl.fr/BALLOON/RUMBA_L2, last access: 5 August 2020. a
Hertzog, A., Vial, F., Mechoso, C. R., Basdevant, C., and Cocquerez, P.: QuasiLagrangian measurements in the lower stratosphere reveal an energy peak associated with nearinertial waves, Geophys. Res. Lett., 29, 70–1–70–4, https://doi.org/10.1029/2001GL014083, 2002. a, b, c, d, e
Hertzog, A., Alexander, M. J., and Plougonven, R.: On the Intermittency of Gravity Wave Momentum Flux in the Stratosphere, J. Atmos. Sci., 69, 3433–3448, https://doi.org/10.1175/JASD1209.1, 2012. a
Hoffmann, L., Hertzog, A., Rößler, T., Stein, O., and Wu, X.: Intercomparison of meteorological analyses and trajectories in the Antarctic lower stratosphere with Concordiasi superpressure balloon observations, Atmos. Chem. Phys., 17, 8045–8061, https://doi.org/10.5194/acp1780452017, 2017. a
Hoffmann, L., Günther, G., Li, D., Stein, O., Wu, X., Griessbach, S., Heng, Y., Konopka, P., Müller, R., Vogel, B., and Wright, J. S.: From ERAInterim to ERA5: the considerable impact of ECMWF's nextgeneration reanalysis on Lagrangian transport simulations, Atmos. Chem. Phys., 19, 3097–3124, https://doi.org/10.5194/acp1930972019, 2019. a
Holt, L. A., Alexander, M. J., Coy, L., Molod, A., Putman, W., and Pawson, S.: Tropical Waves and the QuasiBiennial Oscillation in a 7km Global Climate Simulation, J. Atmos. Sci., 73, 3771–3783, https://doi.org/10.1175/JASD150350.1, 2016. a, b, c
Holt, L. A., Alexander, M. J., Coy, L., Liu, C., Molod, A., Putman, W., and Pawson, S.: An evaluation of gravity waves and gravity wave sources in the Southern Hemisphere in a 7 km global climate simulation, Q. J. Roy. Meteor. Soc., 143, 2481–2495, https://doi.org/10.1002/qj.3101, 2017. a, b
Jensen, E. and Pfister, L.: Transport and freezedrying in the tropical tropopause layer, J. Geophys. Res.Atmos., 109, D02207, https://doi.org/10.1029/2003JD004022, 2004. a, b
Jewtoukoff, V., Hertzog, A., Plougonven, R., de la Cámara, A., and Lott, F.: Comparison of Gravity Waves in the Southern Hemisphere Derived from Balloon Observations and the ECMWF Analyses, J. Atmos. Sci., 72, 3449–3468, https://doi.org/10.1175/JASD140324.1, 2015. a
Kobayashi, S., Ota, Y., Harada, Y., Ebita, A., Moriya, M., Onoda, H., Onogi, K., Kamahori, H., Kobayashi, C., Endo, H., Miyaoka, K., and Takahashi, K.: The JRA55 Reanalysis: General Specifications and Basic Characteristics, J. Meteor. Soc. Jpn. Ser. II, 93, 5–48, https://doi.org/10.2151/jmsj.2015001, 2015. a
Konopka, P., Grooß, J.U., Günther, G., Ploeger, F., Pommrich, R., Müller, R., and Livesey, N.: Annual cycle of ozone at and above the tropical tropopause: observations versus simulations with the Chemical Lagrangian Model of the Stratosphere (CLaMS), Atmos. Chem. Phys., 10, 121–132, https://doi.org/10.5194/acp101212010, 2010. a
Kärcher, B. and Podglajen, A.: A Stochastic Representation of Temperature Fluctuations Induced by Mesoscale Gravity Waves, J. Geophys. Res.Atmos., 124, 11506–11529, https://doi.org/10.1029/2019JD030680, 2019. a, b
Pisso, I., Marécal, V., Legras, B., and Berthet, G.: Sensitivity of ensemble Lagrangian reconstructions to assimilated wind time step resolution, Atmos. Chem. Phys., 10, 3155–3162, https://doi.org/10.5194/acp1031552010, 2010. a
Podglajen, A., Hertzog, A., Plougonven, R., and Žagar, N.: Assessment of the accuracy of (re)analyses in the equatorial lower stratosphere, J. Geophys. Res., 119, 11166–11188, https://doi.org/10.1002/2014JD021849, 2014. a, b, c, d, e
Podglajen, A., Hertzog, A., Plougonven, R., and Legras, B.: Lagrangian temperature and vertical velocity fluctuations due to gravity waves in the lower stratosphere, Geophys. Res. Lett., 43, 3543–3553, https://doi.org/10.1002/2016GL068148, 2016a. a, b, c, d
Podglajen, A., Plougonven, R., Hertzog, A., and Legras, B.: A modelling case study of a largescale cirrus in the tropical tropopause layer, Atmos. Chem. Phys., 16, 3881–3902, https://doi.org/10.5194/acp1638812016, 2016b. a, b, c, d, e, f, g
Podglajen, A., Bui, T. P., DeanDay, J. M., Pfister, L., Jensen, E. J., Alexander, M. J., Hertzog, A., Kärcher, B., Plougonven, R., and Randel, W. J.: Smallscale wind fluctuations in the tropical tropopause layer from aircraft measurements: occurrence, nature and impact on vertical mixing, J. Atmos. Sci., 74, 3847–3869, https://doi.org/10.1175/JASD170010.1, 2017. a, b
Potter, B. E. and Holton, J. R.: The Role of Monsoon Convection in the Dehydration of the Lower Tropical Stratosphere, J. Atmos. Sci., 52, 1034–1050, https://doi.org/10.1175/15200469(1995)052<1034:TROMCI>2.0.CO;2, 1995. a
Preusse, P., Ern, M., Bechtold, P., Eckermann, S. D., Kalisch, S., Trinh, Q. T., and Riese, M.: Characteristics of gravity waves resolved by ECMWF, Atmos. Chem. Phys., 14, 10483–10508, https://doi.org/10.5194/acp14104832014, 2014. a, b
Schoeberl, M. R., Jensen, E. J., and Woods, S.: Gravity waves amplify upper tropospheric dehydration by clouds, EarthSpace Sci., 2, 485–500, https://doi.org/10.1002/2015EA000127, 2015. a
Skamarock, W. C., Snyder, C., Klemp, J. B., and Park, S.H.: Vertical Resolution Requirements in Atmospheric Simulation, Mon. Weather Rev., 147, 2641–2656, https://doi.org/10.1175/MWRD190043.1, 2019. a, b, c
Spichtinger, P. and Krämer, M.: Tropical tropopause ice clouds: a dynamic approach to the mystery of low crystal numbers, Atmos. Chem. Phys., 13, 9801–9818, https://doi.org/10.5194/acp1398012013, 2013. a
Stephan, C. C., Strube, C., Klocke, D., Ern, M., Hoffmann, L., Preusse, P., and Schmidt, H.: Gravity Waves in Global HighResolution Simulations With Explicit and Parameterized Convection, J. Geophys. Res.Atmos., 124, 4446–4459, https://doi.org/10.1029/2018JD030073, 2019. a
Stohl, A., Wotawa, G., Seibert, P., and KrompKolb, H.: Interpolation Errors in Wind Fields as a Function of Spatial and Temporal Resolution and Their Impact on Different Types of Kinematic Trajectories, J. Appl. Meteorol., 34, 2149–2165, https://doi.org/10.1175/15200450(1995)034<2149:IEIWFA>2.0.CO;2, 1995. a, b
Tzella, A. and Legras, B.: A Lagrangian view of convective sources for transport of air across the Tropical Tropopause Layer: distribution, times and the radiative influence of clouds, Atmos. Chem. Phys., 11, 12517–12534, https://doi.org/10.5194/acp11125172011, 2011. a
Ueyama, R., Jensen, E. J., Pfister, L., and Kim, J.E.: Dynamical, convective, and microphysical control on wintertime distributions of water vapor and clouds in the tropical tropopause layer, J. Geophys. Res.Atmos., 120, 10483–10500, https://doi.org/10.1002/2015JD023318, 2015. a
Vincent, R. A. and Hertzog, A.: The response of superpressure balloons to gravity wave motions, Atmos. Meas. Tech., 7, 1043–1055, https://doi.org/10.5194/amt710432014, 2014. a, b
Waite, M. L. and Snyder, C.: The Mesoscale Kinetic Energy Spectrum of a Baroclinic Life Cycle, J. Atmos. Sci., 66, 883–901, https://doi.org/10.1175/2008JAS2829.1, 2009. a
 Abstract
 Introduction
 Datasets and methodology
 Results: intrinsic frequency spectra of Lagrangian fluctuations
 Discussion
 Conclusions
 Appendix A: Aliasing in the frequency spectrum due to the finite output frequency of the reanalyses
 Appendix B: Comparison of GW Lagrangian spectrum estimates using isopycnic, isentropic, kinematic and diabatic trajectories
 Data availability
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Datasets and methodology
 Results: intrinsic frequency spectra of Lagrangian fluctuations
 Discussion
 Conclusions
 Appendix A: Aliasing in the frequency spectrum due to the finite output frequency of the reanalyses
 Appendix B: Comparison of GW Lagrangian spectrum estimates using isopycnic, isentropic, kinematic and diabatic trajectories
 Data availability
 Author contributions
 Competing interests
 Special issue statement
 Acknowledgements
 Financial support
 Review statement
 References