A new approach to global gravity wave momentum flux determination from GPS radio occultatio data

Introduction Conclusions References


Introduction
The important role of gravity waves (GWs) for atmospheric circulation and energy transport is already recognized (Nappo, 2002;Fritts and Alexander, 2003;Sutherland, 2010;Alexander et al., 2010).Recent developments in gravity wave effects and momentum flux analysis are summarized by Alexander et al. (2010).The variety of atmospheric sounding instruments, designed to detect GWs and find answers to related queries, has increased over the last decades.Radiosonde soundings, ground based measurements as well as air-and spaceborne measurements cover a large range of measuring possibilities.The two major measuring techniques Introduction

Conclusions References
Tables Figures

Back Close
Full for satellite measurements are: (i) nadir-scanning instruments, like the Advanced Microwave Sounding Unit (AMSU) (Eckermann et al., 2006) or the Atmospheric In-fraRed Sounder (AIRS) (Hoffmann and Alexander, 2009) and (ii) limb-sounding instruments, like the CRyogenic Infrared Spectrometers and Telescopes for the Atmosphere (CRISTA) (Ern et al., 2004), the High Resolution Dynamics Limb Sounder (HIRDLS), the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) and the many radio occultation (RO) receivers aboard the different Low Earth Orbiting (LEO) satellites like the Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) (see Wright et al., 2011).The GPS RO technique provides vertical profiles of atmospheric properties such as refractivity, from which density (ρ), pressure (P ), temperature (T ) and water vapor pressure (e) are derived.The vertical resolution of RO profiles ranges from 0.1 km in the lower troposphere to 1.4 km in the stratosphere (Kursinski et al., 1997).This technique operates under any weather conditions and provides a global coverage.However, RO events occur without a certain spacing (footprint), neither in time nor in spatial resolution.A long term series of GPS RO measurements is available from the CHAllenging Minisatellite Payload (CHAMP) satellite, which had provided approximately 150 daily profiles between 2001and 2008(Wickert et al., 2001;;Wickert, 2002).A higher density of data is given by COSMIC.The joint Taiwan/US science mission started in April 2006 and delivers approximately 2000 daily profiles (Anthes et al., 2008).RO receivers are also operating aboard of other satellite missions like the Gravity Recovery and Climate Experiment (GRACE, launched in 2002), or the two radar satellites TerraSAR-X (launched in 2007) and TanDEM-X (launched in 2010).The RO technique has been widely used for GW analysis (see Kuo et al., 2000;de la Torre and Alexander, 2005;Alexander et al., 2008a,b;Horinouchi and Tsuda, 2009;Lin et al., 2010;Wang and Alexander, 2010).As shown by Marquardt and Healy (2005), small scale fluctuations of dry temperature RO profiles can be interpreted as gravity waves (GWs), when the vertical wavelength is equal or greater than 2 km.For lower wavelengths the fluctuations may be related to noise.To keep the signal to noise ratio for the temperature fluctuation above the detection threshold, the Introduction

Conclusions References
Tables Figures

Back Close
Full analysis should be applied to altitudes lower than 30 km (Alexander et al., 2008a) or, as suggested by Marquardt and Healy (2005), in the 20-30 km altitude range.The mean specific potential energy (E p ) is a parameter for the characterization of gravity wave activity, given by where N z is the Brunt-V äis äl ä frequency, g is the gravitational acceleration, c p the isobaric heating capacity and T and T are the fluctuation and background components of temperature, respectively.Unlike planetary waves, small scale GWs are usually sub-grid scale phenomena or under-resolved in both the horizontal and the vertical in numerical weather prediction models (Alexander et al., 2010).In order to quantify the effect of small scale GWs alone, it is important to separate these from planetary waves.The GW analysis is based on the extraction of the vertical T contribution.Together with the measuring technique, the background determination sets the range of detectable GWs even though the amount of planetary waves can not be perfectly extracted.Here, a continuous wavelet transform (CWT) analysis is used for the background determination (see Sect.2.1).Ern et al. (2004) introduced a method to derive horizontal wavelength (λ p h ) projected along the line connecting two adjacent T profiles (using CRISTA temperature data), for a given altitude: where ∆x i j is the distance and ∆φ i j is the phase shift between the two profiles at a given altitude, respectively.To extract λ h , at least a third T profile is needed (not Introduction

Conclusions References
Tables Figures

Back Close
Full aligned with the others).Measurements must be close in time and space, to provide observations of the same wave packet.The absolute value of vertical flux of horizontal momentum (MF), for mid frequency approximation and assuming a single prevailing wave within the packet, is given by where λ v is the vertical wavelength.Wang and Alexander (2010) used clusters of three or more RO measurements grouped in 15

GPS RO technique
A GPS receiver on board a LEO satellite measures the Doppler shift of the dual frequency signal of one of the many GPS satellites rising or setting on the local horizon.The Doppler shift can be converted to the radio wave bending angles from which the atmospheric refractivity is derived.In the neutral atmosphere, the refractivity is further reduced to temperature (T ), pressure (P ) and electron density profiles (Anthes et al., 2008;Kursinski et al., 1997).Introduction

Conclusions References
Tables Figures

Back Close
Full The capability of the RO technique to properly detect GWs depends principally on the relative orientation between the wave fronts and the line-of-sight (LOS) between the GPS and LEO satellites.Wave amplitudes can be better resolved when the fronts are nearly horizontal or when the angle between the LOS and the horizontal component of the wave vector approaches π/2 (Alexander et al., 2008a).Short horizontal scale waves have a high probability to be attenuated or not be detected at all.In opposite to CRISTA, the GPS RO data sets (indifferent of which mission) have neither a regular spanned grid on which the measurements are taken nor an almost constant LOS.
In this work we use the COSMIC pos-processed dry temperature profiles provided by the UCAR COSMIC Data Analysis and Archival Center (COAAC) Version 2007.3200 for the global analysis.The dry temperature profiles are provided from near the surface up to 40 km, at 0.1 km steps.

Data processing
In order to calculate E p , the fluctuation has to be extracted from each T profile.Earlier works (Hocke and Tsuda, 2001;Tsuda et al., 2000Tsuda et al., , 2004;;Baumgaertner and McDonald, 2007;de la Torre et al., 2006a,b;Schmidt et al., 2008) used a vertical filter directly on the T profile.The resulting T profiles (T = T − T ) still include longer waves with vertical wavelengths (λ v ) similar to GWs (e.g.Kelvin waves Holton et al., 2001).To exclude planetary and Kelvin wave contribution, a T including these waves must be computed.Following closely Wang and Alexander (2009Alexander ( , 2010) ) we first binned all T profile data (taking into account a compromise between the resolution that we may obtain and the density of available profiles) into a 10 × 15 • lat/lon field between 9-40 km with a vertical resolution of 100 m.For the daily binned profiles, a CWT as a function of longitude is performed, for each latitude and altitude.This allows to estimate wave numbers 0-6 for each longitude.The "large scale variation", reconstructed by wave numbers 0-6 are then interpolated back onto the location of the original T profile and subtracted from it.The resulting T represents the wave fluctuations of short scale waves (shorter than Introduction

Conclusions References
Tables Figures

Back Close
Full wave number 6).For this determination of T time variations of large-scale waves within 1 day are not taken into account.
Then λ v (using a CWT analysis) and E p (using Eq. 1) can be directly computed from the T profiles and are later presented as seasonal mean values in a 5 × 5 • resolution.From Eq. ( 2), λ p h along the line connecting two adjacent T profiles for a given altitude is determined.The method is restricted to pairs of measurements within a small temporal and spatial window (both measurements are supposed to be taken within the same π-cycle of the wave).To ensure the same wave packet, the vertical structure of the considered profiles must be similar.For the phase difference determination, it is very important that all the RO profiles are measured with a similar LOS.In ideal case, the LOS of all the measurements are parallel.In addition to that, the geometric configuration between the LOS and the wave propagation direction limits the capability to see a GW (Alexander et al., 2008a).Using at least three RO measurements (not collinear with each others), we can derive λ h instead of λ p h .All RO profiles of these groupings must fulfill strict temporal and spatial restrictions.In this work the temporal and spatial windows are set to 2 h and 15 • in all directions.The reason for these parameter settings will be explained in Sect. 4.
Alexander and de la Torre (2010) used a similar method to derive wavelengths in a Cartesian coordinate system from three soundings.They describe that the ideal case would be a perpendicular arrangement between the three soundings, but the COSMIC constellation does not provide these ideal cases.
The six identical COSMIC satellites were launched on one rocket into an orbit of 405 km and later on raised to their final orbit of approximately 800 km.During the COS-MIC clustering phase (the first months) several measurements were obtained in a small spatial and temporal window and with similar LOS.
Before this method is applied to the RO data, the relative positioning of pairs of RO profiles with respect to the wave must be taken into account.In the case that the alignment of one of the pairs is perpendicular to the wave propagation direction, the T profiles will not have a phase shift between each other.Those pairs of RO measurements, Introduction

Conclusions References
Tables Figures

Back Close
Full where ∆Φ is lower than 0.5 radian (λ p h > 10 000 km), are not taken into account (see Sect. 4).
From the CWT analysis (using a Morlet wavelet) the dominant λ v at considered altitudes is determined.We treat T profiles with a dominant λ v within an interval of 2 km as they belong to the same wave system.Ern et al. (2004) used a window of 6 km for the CRISTA-2 data set.In Fig. 1a, two RO temperature profiles (T 1,2 ) are displayed for the altitude range between 20-30 km.Their global (for the entire altitude range) dominant λ v is shown in Fig. 1b for T 1 and 1c for T 2 .The cross-wavelet spectrum (Fig. 1d) shows the amplitudes of the dominant λ v at each altitude.This is needed for the phase shift determination at each altitude between all pairs of T profiles.Once the phase shift and the distances are known, k p h (k h = 2π/λ h ) can be derived for all altitudes.

Horizontal wave parameters
For the MF determination, we need k h instead of k p h .As mentioned before, a third measurement is needed to derive k h .For a detailed explanation we take a look at Fig. 2, starting with 2a.Two vertical T RO profiles (T 1 in brown and T 2 in blue), separated by the horizontal distance ∆x are displayed.The ∆Φ between these measurements is derived from a cross wavelet analysis for each altitude.Figure 2b, shows the derived k p h (orange arrow) from three T measurements at one altitude and the true k h (green arrow).The dotted lines are the corresponding λ p h (orange line) and λ h (green line).When using these three measurements, k h can be derived by: where (k, l ) are the horizontal components of the wave number vector, and (x i ,j , y i ,j ) are the RO coordinates at given altitude.∆Φ i j is derived from cross wavelet analysis for each pair of T profiles (Fig. 2c).Since this is an over-determined problem, there can Introduction

Conclusions References
Tables Figures

Back Close
Full be more than one result.The reference point (RP) from which this system of equations is solved, plays an important role.As shown in Fig. 2c, the location of (x i j , y i j ) and the distances between the points (∆x i j ) are known, therefore we can derive the phase shifts (∆Φ i j ) between each measurement pair.When using P1 as RP, then Eq. ( 4) is solved using two equations for ∆x i j : ∆x 12 and ∆x 13 .In a similar way, P2 as RP ∆x i j are ∆x 12 and ∆x 23 and for P3 as RP ∆x i j are ∆x 13 and ∆x 23 .When solving the equations for each RP, three k h are derived.As displayed in Fig. 2d at least two identical k h (green arrows) and a third k h (purple arrow) are derived.All three displayed green arrows are identical in lengths and orientations.The two starting in P1 and P3 represent the computed wave vectors, the third green arrow is shifted to fit the phase maximums of the background wave.The inconsistent k h can be different from the two identical k h in both, orientation as well as length of k h , but it might also be identical to the two others.This depends on the angles between the three measuring points as explained below.When taking a look at the wave fronts the reason for the wrong determination of k h at one RP becomes more clear.The RP that produces the non consistent result is placed in between the other two points regarding the wave fronts.
Another way to treat this problem is a geometric description as presented in Figs. 3  and 4. Considering the same three points as before, representing the intersection of three closely separated RO profiles with a constant altitude plane.The unknown constant successive wave phase lines Φ and Φ ± 2π are illustratively shown in Fig. 3a.
The line of action (line ab), as well as the direction of k, are obviously also unknown.Now, the three possible pairs of points, as shown in Fig. 3b-d are considered.The slope of k h may be determined, for example in the case shown in Fig. 3b, in the following way: from Eq. ( 2) and the knowledge of both the phase and distance differences between each pair of points, λ p h (black sinusoidal lines) may be determined (Fig. 4a).If now we displace these curves along the blue lines defined by each pair, until matching a common extreme (red sinusoidal curves), the opposite extremes of these red curves define the slope of Φ (line cd).We repeat this procedure in Fig. 4b, observing identical result.In the case of (Fig. 4c), note that one of the black curves was displaced in the Introduction

Conclusions References
Tables Figures

Back Close
Full opposite direction, in order to be able to obtain the same slope as in the precedent two cases.When displacing the black line connection the brown and blue points as before (gray sinusoidal line) the resulting slope is presented by the gray line between ef, which is not confirm with the results of Fig. 4a, b.This happens, when one of the inner angles α, β or γ (in Fig. 3a) is an obtuse angle.If the three angles are acute, no inversion is needed and the three possibilities described, render straightforwardly an identical result.
In general, for practical purposes and to handle a huge amount of triads of points, after computing the phase slope as in Fig. 4a-c for each triad, no possible uncertainty remains if we simply consider the prevailing slope (to include triads with one inner obtuse angle).λ h is finally illustrated in Fig. 4a-c (green curves).

Sensitivity study
As mentioned in Sect. 1 the temporal window should be kept as small as possible, still allowing a sufficient number of triads to extract global distributions of the wave parameters.We choose a 2 h interval.Shorter windows would decrease the number of possible triads too much, whereas a larger time window would increase the number of triads.But ∆x i j will be unrealistic because of the phase shift variations from profiles from different wave origin would be included.For a closer look, we apply a sensitivity study to Eq. ( 2), using two sinusoidal curves with different phase shifts (0 • and 180 • ) and four horizontal distances (100, 200, 500 and 1000 km).This only focuses on the one-dimensional case (on λ p h ).With a decreasing phase shift, λ h approaches very high values (Fig. 5).This effect is even stronger when ∆x increases.A detection of small phase shifts (lower than 30 • ) shows for all distances very long wavelengths (λ p h → ∞).For larger distances (for example 1000 km), a phase shift of 45 • corresponds to a λ p h of more than 10 000 km.These long λ h are unrealistic for GWs.This and the fact, that alignments cannot be ruled out (small phase shift at a given altitude might be due to the ±2π or even ±3π periodicity), leading to a limitation of the phase shift.Therefore we Introduction

Conclusions References
Tables Figures

Back Close
Full set a lower limit of 0.5 radian (∼ 30 • ), when additionally the maximum distance is set to be 15 • .The new limiting factor must be taken into account when using this method to analyze GPS RO profiles.Now we take a look on the impact of the distance dependence for the λ h determination using Eq. ( 4) (the two-dimensional case).Triads with the spatial interval of 5 • , 10 • and 15 • were identified.The first and obvious result is that the shorter the interval, the less triads can be found.This was performed for the Northern Hemisphere (NH) winter 2006/07 for the altitude range 20-30 km.For 5 • the amount of triads is too small for a sufficient statistical global analysis.The results for 10 and 15 • are shown in Fig. 6.
Although the absolute values vary by approximately 50 %, the structure of λ h is quite similar for both settings.The same effect can be found when studying the results for the MF distribution.A more detailed description of λ h and MF is given in Sect. 5.This method can provide λ h only as an upper boundary due to alignment effects and limitations in the measuring principles and data availability.The 15 • horizontal spacing for the triad search was chosen for the global analysis.In addition, the data availability does not allow a global analysis for NH summer ( 2006) with a 10 • interval.(see Fig. 7c, d) is similar, but the clustering phase of the beginning of the COSMIC mission has a small impact which leads to shorter distances for the NH summer than winter.These differences are most distinct over Australia and Southern America.The maximum distance is limited by a 15

Global analysis
• circle, which represents about 1500 km at the equator and 500 km at high latitudes.The gaps in the NH summer are found in all wave parameters.There, the number of triads is to small for a statistical analysis.For all wave parameters, as well as for the mean distance, the results of each triad (mean E p , λ v and the computed λ h along with MF) are presented for their mean lat./lon.position (MP) and then fed into a 5 × 5 • lat./lon.grid.Due to the low numbers of triads in polar regions (north of 75 • N, south of 75 • S) these areas are not evaluated.As mentioned before, a temporal and spatial window is used, the ∆Φ limitation and the comparison of the dominant λ v eliminate about 30 % of all measured triads.Figure 8a, b displays the λ v distributions for NH summer (a) and winter (b) with values of minimum 5 km around the equator and up to 10 km in the extra tropics and polewards.At first sight, a small general decrease of λ v from Fig. 8a to b can be seen.
Taking a closer look, some significant changes in both directions are detectable.For example the decrease of λ v over Canada and the US as well as over Russia from about 9 km down to less then 7 km or the increase of λ v over Australia and north Africa.All these changes are located over land, but there is also a decrease in the λ v over the Southern Hemisphere (SH) water regions, which is most significant around the Antarctic Peninsula.
The E p distribution (see Fig. 8c, d) agree with previous results (Schmidt et al., 2008).the north east of Russia, the NH winter weather conditions provide strong wind shear.
In northern Europe, the prevailing westerlies induces GWs on the lee side of the Scandinavian mountains with values of up to 5 J kg −1 .The subtropical jet stream induces hight values of E p when crossing the Zagros Mountains (Farajzadeh et al., 2011).For both seasons, the equator region shows high E p values around 4 J kg −1 .These values are due to the high amount of convective potential from high temperatures and humidity in that area.
The λ h distributions, see Fig. 9a for JJA and 9b for DJF, show a trend of increasing λ h towards the equator.The reason for the different starting point of the λ h scale is the visibility of the pattern which fads in other scales for each graphic.This kind of pattern has already been shown by Fr öhlich et al. (2007), who used a theoretical λ h , and by Wang and Alexander (2010) using GPS RO data.Also, the λ h distribution shows longer λ h over water on the winter hemispheres than on the summer hemisphere and decreasing λ h over the winter hemisphere land areas compared to the same regions in summer.For the SH winter (Fig. 9a) the effect of shorter wavelength over land is most significant.The absolute value of λ h at the east of the Andes is higher than the results from de la Torre and Alexander (2005), but the here derived λ h is only an upper boundary of λ h (see Sect. 4).The λ h would decease by a mean of 50 % when changing the spacial window from 15 to 10 • .
The global MF distribution is displayed in Fig. 9a, b for JJA and DJF respectively.The absolute values of MF range between 0.2 and 1.7 mPa, but the pattern of the distribution is more crucial for the evaluation.The absolute values depend on the derived λ h and therefore also in the chosen spacial interval.For the MF however, the absolute values increase when choosing a shorter spacial window, but the pattern does not change significantly.For both seasonal mean MF distributions (JJA and DJF), a band of high values is located along the equator.Additional to that, the winter hemisphere shows more wave activity (higher MF values) than the summer hemispheres.The higher values are found at the east of the Andes mountains and along the Antarctic Peninsula Introduction

Conclusions References
Tables Figures

Back Close
Full during the SH winter (Fig. 9c) and at the east on the Scandinavian mountains and the Zargos mountains (Iran) during the NH winter (Fig. 9d).
Other global distributions of the MF were given, e.g. by Ern et al. (2004Ern et al. ( , 2011););Fr öhlich et al. (2007); Wang and Alexander (2010).Again, the absolute values shall not be compared even though the range is in the 0.2 to 1.2 mPa fits the results derived by Ern et al. (2004) andFr öhlich et al. (2007).But the general patterns fit each other as well.Not all results are presented for both seasons and for that exact time period, therefore the patterns can also differ.But all results show higher values of MF on the winter hemisphere.Wang and Alexander (2010) show also the two spots over southern America and Africa for the same time period (DFJ) and the high MF over Scandinavia.
The method of taking triads instead of a fixed 15 × 15 • grid, allows a higher resolution in the final distributions of the derived wave parameters.

Conclusions
This work proves that the derivation of MF from GPS RO data is possible when certain restrictions and satellite constellations are regarded.COSMIC data from the clustering phase (April 2006-February 2007) were analyzed.Some of the limitations applied come from the method introduced by Ern et al. ( 2004) (small temporal and spatial windows, analysis for each altitude), others were added to minimize the alignment effects ( λ v , ∆Φ and E p ). Global distributions of the vertical, horizontal wavelength as well as potential energy and MF as seasonal means in the altitude range 20 to 30 km, from June-August 2006 and December 2006-February 2007 are displayed as an example of our studies.The new applied phase shift limitation and the restrictions in the geometric constellation of the regarded triads lead to a higher resolution in the results then Wang and Alexander (2010) showed.However, a disadvantage of this method is the high dependency of λ h on the spacing within the triads.After the clustering phase of the six COSMIC satellites, the occultations were further apart making it harder to find Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

For a global
analysis, seasonal mean data for June-August 2006 (JJA) and December 2006-February 2007 (DJF) are shown.The altitude range for the statistical analysis is set to 20-30 km.The global distribution of the number of available triads and their mean distance between each triad are shown in Fig. 7 for the NH summer (a, c) and winter (b, d), respectively.Figure 7a, b show, that the number of triads are higher in the extra tropics then in the equator region or closer to the poles.As mentioned in Sect.4, the amount of measurements in the beginning of the COSMIC mission was not as high as expected.By the end of 2006 all six satellites were running and therefore increasing the number of ROs.For Fig. 7a, b, different scales are chosen to see the structure of the grouping of triads for both seasons.The mean distance of the triads Discussion Paper | Discussion Paper | Discussion Paper | The higher values (up to 5 J kg −1 ) are located along the equator for both seasons and decrease towards the polar regions.However for June-August 2006 the higher values in the extra tropics are located on the SH, with explicit high values along the east side of the Andes mountains mainly due to low level westerlies and the permanent jet (Alexander and de la Torre, 2010).For December 2006-February 2007 higher values are found in the NH.Going from west to east, the increase in E p in the NH can be explained by the various mountain chains and wind shear conditions.East of northern America, and in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | suitable triads.But, the method can be applied to other vertical temperature measurements, as long as they provide the required spatial and temporal resolution.Discussion Paper | Discussion Paper | Discussion Paper | Hocke, K. and Tsuda, T.: Gravity waves and ionospheric irregularities over tropical convection zones observed by GPS/MET radio occultation, Geophys.Res.Lett., 28, 2815-2818, doi:10.1029/2001GL013076,2001.2912 Hoffmann, L. and Alexander, M. J.: Retrieval of stratospheric temperatures from Atmospheric Infrared Sounder radiance measurements for gravity wave studies, J. Geophys.Res., 114Discussion Paper | Discussion Paper | Discussion Paper | Tsuda, T., Nishida, M., Rocken, C., and Ware, R. H.: A global morphology of gravity wave activity in the stratosphere revealed by the GPS occultation data (GPS/MET), J. Geophys.Res., 105, 7257-7273, doi:10.1029/1999JD901005,2000.2912 Tsuda, T., Ratnam, M. V., May, P. T., Alexander, M. J., Vincent, R. A., and MacKinnon, A.:Characteristics of gravity waves with short vertical wavelengths observed with radiosonde and GPS occultation during DAWEX (Darwin Area Wave Experiment), J. Geophys.