Temperature Profiles From Two Close Lidars and a Satellite to Infer the Structure of a Dominant Gravity Wave

Gravity waves (GW) are a crucial coupling mechanism for the exchange of energy and momentum flux (MF) between the lower, middle, and upper layers of the atmosphere. Among the remote instruments used to study them, there has been a continuous increment in the last years in the installation and use of lidars (light detection and ranging) all over the globe. Two of them, which are only night operating, are located in Río Gallegos (−69.3° W, −51.6° S) and Río Grande (−67.8° W, −53.8° S), in the neighborhood of the austral tip of South America. This is a well‐known GW hot spot from late autumn to early spring. Neither the source for this intense activity nor the extent of its effects have been yet fully elucidated. In the last years, different methods that combine diverse retrieval techniques have been presented in order to describe the three‐dimensional (3‐D) structure of observed GW, their propagation direction, their energy, and the MF that they carry. Assuming the presence of a dominant GW in the covered region, we develop here a technique that uses the temperature profiles from two simultaneously working close lidars to infer the vertical wavelength, ground‐based frequency, and horizontal wavelength along the direction joining both instruments. If in addition within the time and spatial frame of both lidars there is also a retrieval from a satellite like SABER (Sounding of the Atmosphere using Broadband Emission Radiometry), then we show that it is possible to infer also the second horizontal wavelength and therefore reproduce the full 3‐D GW structure. Our method becomes verified with an example that includes tests that corroborate that both lidars and the satellite are sampling the same GW. The improvement of the Río Gallegos lidar performance could lead in the future to the observation of a wealth of cases during the GW high season. Between 8 and 14 hr (depending on the month) of continuous nighttime data could be obtained in the stratosphere and mesosphere in simultaneous soundings from both ground‐based lidars.


Introduction
Gravity waves (GWs) have significant global effects from the lower to the upper atmosphere (e.g., Fritts & Alexander, 2003;Gill, 1982). They are mainly generated in the troposphere or stratosphere and may increase in amplitude under vertical propagation in certain conditions. These waves may transfer significant amounts of energy and momentum flux (MF) to the background if filtering or dissipation occurs while they propagate. This may result in strong forcing of the dynamics and thermal structure, mainly in the middle atmosphere. Some works have shown that GW may even penetrate and influence the thermosphere and ionosphere (e.g., Alexander et al., 2015;Park et al., 2014).
Although the austral tip of South America may be the most intense hot spot of GW on the globe from austral late autumn to early spring (e.g., Ern et al., 2004;Hoffmann et al., 2013Hoffmann et al., , 2016, several studies using lidars (light detection and ranging), aircrafts, radars, or balloons have focused on the Northern Hemisphere (NH). However, in the last years there has been a growing awareness of the relevance of this GW hot spot close to the southern pole (e.g., Chu et al., 2018;Fritts et al., 2016;Kaifler et al., 2015;Llamedo et al., 2019;Zhao et al., 2017). The importance of an improvement in the knowledge of this region is highlighted by the fact that comparisons of stratospheric GW MF obtained from general circulation models (GCMs) and satellite data reveal some notable discrepancies. Although there are some well reproduced features, large deviations are still present (e.g., de la Cámara et al., 2016;Geller et al., 2013). For example, several GCMs produce simulations of the Southern Hemisphere (SH) polar stratosphere that lead to significant underestimations of the temperatures and drag on the winds (e.g., Butchart et al., 2011;Wright & Hindley, 2018). Numerical solutions are not able to resolve the full spectrum of waves. Parameterizations of the smallest scale GW are then introduced, but they are usually too coarse and may be a major cause of the biases in the polar SH stratospheric dynamics and thermal structure simulations (McLandress et al., 2012). These shortages recall the need for observational information on GW sources, evolution, and general behavior in this zone. The possible but still uncertain causes of the intense GW activity are usually attributed to orography (Southern Andes, Antarctic Peninsula, or small oceanic islands), nonorographic waves from winter storm tracks over the Southern oceans or from spontaneous adjustment, or jet instability around the edge of the stratospheric vortex or secondary waves stemming from primary breaking ones from any source (e.g., Alexander & Grimsdell, 2013;Hindley et al., 2015;Sato et al., 2009).
The DIAL (Differential Absorption Lidar) instrument, belonging to the Centro de Investigaciones en Láseres y Aplicaciones (CEILAP), was located in 2005 at the Observatorio Atmosférico de la Patagonia Austral (OAPA) in Río Gallegos (51.6 • S, 69.3 • W), mainly for ozone studies. This lidar was the southernmost to the North of Antarctica until November 2017, when CORAL (Compact Rayleigh Autonomous Lidar) started working in Río Grande (53.8 • S, 67.8 • W). Río Gallegos is nearly 300 km to the East of the Andes and 70 km to the North of the Strait of Magellan, whereas Río Grande is further Southeast on the Atlantic coast of the Tierra del Fuego island (Figure 1). Both are in an excellent position in relation to the observation of the GW hot spot and are separated by 265.6 km (zonal and meridional distances of 100.8 and 245.7 km).
It is not possible to determine the MF of a dominant GW and its 3-D structure with the temperature retrieval from one lidar. However, it may be feasible to obtain additional information with a second simultaneous and close lidar through the phase shift between GW-induced perturbations on both soundings and the knowledge of the spatial separation of both instruments. If in addition there is another close temperature profile, for example, provided by a satellite, then it may be possible to reveal the full 3-D GW structure, including the net MF calculation. The satellite measurements like those from SABER (Sounding of the Atmosphere using Broadband Emission Radiometry), GPS radio occultation, or HIRDLS (High Resolution Dynamics Limb Sounder) provide no directional information of the horizontal components of the MF vector, but only the absolute value of each one can be found with 3 near profiles. However, an inspection of the horizontal 10.1029/2020EA001074 components of the equation of motion for the atmosphere shows that the net MFs affect the wind and temperature structure (e.g., Geller et al., 2013). Wright et al. (2016) also remarked the importance of obtaining the net rather than the absolute value. In brief, we suggest to use a sequence of vertical temperature profile pairs over a time interval plus a static retrieval. As far as we know, there have been no previous similar studies. Just frozen GW reconstructions have been usually obtained from a combination of instantaneous satellite temperature profiles, which are close in space and time (e.g., Alexander et al., 2018;Ern et al., 2004Ern et al., , 2017Schmidt et al., 2016), as the evolution could not be monitored in those cases.
In the present study we employ a succession of two close and simultaneous lidar temperature measurements over time and height and a third instantaneous retrieval from a satellite within the same spatial and time frame in order to infer the ground-based frequency and the three cartesian wavelengths of a dominant GW in the studied zone. This also allows the determination of the three phase velocity components and the net GW MF. Section 2 explains the analysis procedure here employed and the general characteristics of the data from both lidars and from the SABER instrument onboard the TIMED (Thermosphere Ionosphere Mesosphere Energetics Dynamics) satellite that we use in an application example in section 3. Section 4 summarizes the constraints of the method and the main results of our case study.

Data and Method
The DIAL lidar at Río Gallegos has four Newtonian telescopes of 0.5 m diameter. Four Rayleigh and two Raman digital channels record the backscattered photons emitted by the third harmonic of the Nd-YAG laser at 355 nm (130 mJ maximum energy) with a 30 Hz repetition rate. For details on DIAL working characteristics, see Llamedo et al. (2019). Above 30 km altitude, temperature T is obtained using the Rayleigh scattering technique, whereas below the method is affected by aerosol scattering and ozone absorption. The spatial/temporal resolution of the photon counting system is 15 m/1 min, respectively, but an integration of at least 900 m/30 min is needed to improve the signal-to-noise ratio (SNR). A careful analysis revealed that a significant fraction of oscillations above 40 km height with a 30 min integration time may be caused by noise. Moreover, very large negative temperature gradients are unlikely to persist as they are convectively unstable. The DIAL signal power is about 4 W.
The CORAL lidar measures atmospheric backscatter profiles from 22 to 90 km altitude, but only values above 30 km should be used due to the effect of aerosols. It is an Nd:YAG laser generating 12 W at 532 nm and 100 Hz pulse repetition rate. The telescope comprises a 630 mm diameter f/2.45 parabolic mirror. At the top altitude, the T derivation procedure is seeded by SABER temperature. Retrievals with a temporal and vertical resolution of, respectively, up to 10 min and 0.3 km may be provided with reasonable SNR values. For a description of CORAL characteristics see Kaifler et al. (2017).
With the two lidars we may obtain a two-dimensional (2-D) scenario over several hours. In order to be able to fully resolve the 3-D structure of a GW observed by both lidars, an additional profile must be provided. SABER soundings (Mlynczak, 1997), if present within the time and space frame, seem to be a good option as they measure T roughly between 20 and 100 km height. Although this retrieval is instantaneous, it helps to resolve the 3-D GW structure over the whole observational period of both lidars (assuming the wave persists and does not suffer a substantial modification during all that time). Here we use kinetic temperatures from Version 2.0 datasets. A first step requires the separation in all the T profiles of GW from the background, including planetary waves (PWs) if present. In general, special care has to be taken in avoiding spurious amplifications of GW near the altitudes of sharp changes (tropopause or stratopause) when using a digital filter to isolate these waves (e.g., de la Torre et al., 2006). Due to the limited vertical range of reliable lidar data in Río Gallegos (30-40 km), we do not undergo that problem. We follow Ehard et al. (2015) and Rapp et al. (2018) in that at middle and high latitudes a digital filter cutoff at 15 km of vertical distance separates GW from PW. As we have a 10 km vertical interval of data, a band pass between 2 and 10 km will, respectively, contemplate the Nyquist condition (a vertical resolution of 1 km and sampling of 100 m is used for consistency in both lidars) and the elimination of the background and PW. We used a Savitzky-Golay filter (Orfanidis, 1996). Regarding the time evolution of both lidars, we implemented for consistency intervals of 30 min. The same filter was used below for the SABER profile.

10.1029/2020EA001074
We assume that a dominant monochromatic GW is present at both lidars during the whole sounding period or at least a significant fraction of it. In the last case the different portion should be identified from the data. The perturbed temperature T' A,R , respectively, at Río Gallegos ( A ) and Río Grande ( R ) may be then represented by where x, y, z, and t represent zonal, meridional, and vertical coordinates and time; k, l, m, and are the corresponding wave numbers and frequency as seen from the ground; T o is the GW amplitude (we consider it to be constant due to the limited height range); and 0 is a fixed value, whereas on the right-hand side the expression within parentheses is the wave phase. If for both lidars we represent as usual T ′ against z and t, we should then get a similar and m if both places are observing the same dominant wave and are subject to similar mesoscale conditions. If so, it means we obtained the ground-based frequency and the vertical wavelength.
Between the two lidars at a fixed time and height, the phase difference d AR is given by If we put the origin in Río Grande and rotate the horizontal cartesian coordinate system so that the y * axis coincides with the direction to Río Gallegos (see Figure 1), then we may rewrite Equation 3 shows that if the phase difference is found, as the horizontal separation between both lidars is known, it is possible to obtain the component of the horizontal wave vector defined by the direction that joins both places. To reconstruct the full 3-D GW structure, only one horizontal wavelength is missing. This information may be provided by an additional profile. To optimize its added value, its location in the horizontal plane should not lie along the y * axis (the new information would be redundant) but rather separated from it.
In addition, to ensure that it may be observing the same dominant GW, we restrict its separation from any of both places to less than 2.5 • in latitude and 4 • in longitude (the angle difference keeps the maximum possible zonal and meridional separations equal). A minimum distance, for example, 50 km, should be set to avoid uncertainties being larger than the possible small phase difference for a too close comparison. It will be shown below that the SABER horizontal excursion for measurements between 30 and 40 km height is small compared to the horizontal wavelength found in our example, so we essentially consider it a vertical profile. The second horizontal equation between SABER and Río Gallegos then is where l * was already found in Equation 3, the positions are known, so k * can be calculated.
The determination of GW phase differences between both lidars as a function of height at consecutive 30 min intervals has been performed by wavelet coherence (Torrence & Compo, 1998). We first assumed that d AR was between − and . However, we then also contemplated the three other possible aliased phase differences between −4 and 4 . This implies physically that we also contemplate waves with smaller wavelengths and/or propagating in the opposite direction than the initial one. Aliasing is not expected in the vertical direction or time, as it would imply wavelengths smaller than 2 km and ground-based periods of less than 1 hr. We will then obtain four possible solutions. Once we derive the missing horizontal wavelength for each of the four aliased cases, the "true" value will be selected as the one that better suits the GW dispersion relation. It is given in terms of the intrinsic frequencŷbŷ where H is the scale height, k 2 h = k 2 + l 2 = k * 2 + l * 2 is the squared total horizontal wavelength and N is the Brunt-Väisälä frequency, which will be obtained from ERA Interim reanalysis T profiles adequately interpolated in time and space but may be also obtained from any lidar (Chu et al., 2018). The reanalysis 10.1029/2020EA001074 horizontal components of wind were used in order to calculate the intrinsic frequency on the left-hand side (̂= − kU − lV with U and V the zonal and meridional projections of air velocity). A similar procedure to discard spurious aliased cases was already used by Alexander et al. (2018).
Fortunately, the part of the year with the longest nights coincides with the months of the most intense GW. From March to October 2018 there were 17 coincident measurement periods of both lidars, ranging from 4 to 12 hr of simultaneous sounding intervals. The relatively low number of concurrent observations as compared to the total number of nights is due to cloudy conditions or to operational problems in any of both places. To make both data sets comparable, the information from Río Grande was restricted to the 30-40 km height interval, all the profiles were provided every 30 min, and the vertical resolution was set to 1 km. The T retrievals were initially used from 27 to 42 km height, but after retaining the GW, they were restricted to 30-40 km. This procedure was only done to attenuate any artificial discontinuity at the beginning or end of the data set due to the implicit assumption of the digital filtering procedure that it is cyclic, which may introduce spurious temperature fluctuations (Ehard et al., 2015). Although all the profiles already underwent quality control verifications, we tested them against anomalous temperature values (below 160 K or above 320 K).
After the 17 matrix pairs of T ′ against height (every 100 m in the 30-40 km range) and time (every 30 min along the coincident observational period) were obtained, different procedures were developed to ensure that both lidars detect the same GW and adequately. A minimum of 10 consecutive with up to one missing measurement (to be interpolated) time was requested. By visual inspection we kept only cases that exhibited the presence of wavefront-like features and eliminated cases with clearly identifiable noisy patterns in either lidar site, whereby eight cases out of 17 passed this selection process. In every matrix pair we searched at every fixed time a dominant vertical mode that was present and found its phase difference at both places by wavelet coherence (Torrence & Compo, 1998). If d AR had an abrupt change at any height (>0.2 rad in 100 m), then the pair was discarded as it may mean that different phenomena were observed at both places. Or diffusion or absorption or any instability happened or simply the main observed effects cannot be explained in terms of a significant GW or its properties stayed beyond the observational window of the lidars (we recall that no instrument may capture the whole spectrum of waves). With these requirements only four nights from March to June were still suitable for further analysis.
For these remaining cases we applied 2-D wavelet analysis to each T ′ (t, z) matrix (Wang & Lu, 2010;Kaifler et al., 2017). The spectral power (SP) as a function of height, time, vertical wavelength, and period was obtained for each lidar and event, whereby the largest value indicated the dominant mode and location. In order to ensure that the main GW seen at each site was similar to the other one, we required that they both become represented by slightly differing elements of the 2-D wavelets basis: the angle and magnitude that define the dominant mode for each lidar sounding should deviate by less than ∕10 and by less than 10 units (as vertical spacing is 100 m, this represents, e.g., 1 km in the z direction). After this evaluation only two cases were found to meet these criteria. However, when requiring that the evaluation of Equation 5 should differ by less than 10% as calculated on both hand sides, only one event remained. Larger deviations could mean that we are not observing a GW or it may be undergoing nonlinear behavior at any of both places. We show in the appendix some characteristics of the case that missed our last test and give some remarks on its failure.

Application Example
The case to be analyzed is 1 June 2018 21:50 to 2 June 2018 10:08 (all times in UTC). In Figure 2 we show SP for both lidars. The outcome of the 2-D Morlet wavelet is determined in terms of two parameters: the angle , which defines the direction of the mode in t − z space, and the scale s, which is the wavelength along the angle direction. SP was initially obtained as a function of z, t, , and s, whereby the summation over the former two variables yielded SP as a function of the latter two. Differences in angle and scale for the dominant mode in both lidars were, respectively, less than ∕20 and 0.1 km, and their values were ∕2 and 5 km. Notice that = ∕2 implies stationary wavefronts, as they are represented by horizontal lines in the t, z plane. In Figure 3 we show the location of the polar vortex. To obtain it at different heights, we used ERA Interim data at 475, 600, and 700 K isentropic levels to find the largest potential vorticity gradient weighted by the horizontal wind speed (Nash et al., 1996). It can be seen that both lidars are outside and far away from the vortex, so direct effects on both soundings or the presence of nonstationary GW induced by geostrophic adjustment are unlikely. It should be mentioned that in some occasions the edge may reach or even surpass one or both sites. In addition, to analyze if all the studied region exhibits similar mesoscale features, we show in Figure 4 the horizontal velocity at 600, 100, and 10 mb levels. Although the dominant wind direction changes with height, homogeneous conditions can be observed at every single altitude shown. We therefore expect uniform background characteristics in the whole observed zone. In addition, in the lower left part of the 600 mb panel (a level close to the height of the local Andes mountains), adequate conditions for the generation of mountain waves are observed. The prevailing wind finds a meridional obstacle, thus probably generating GW (Baines, 1995). Moreover, the other panels show that at least at those two other heights, critical levels for stationary mountain waves are unlikely (no zero wind regions).
In Figure 5 we show the T ′ representation against time and height for both lidars. Notice some general similarities regarding the nearly stationary wavefronts. The best fit wavefront maxima are also shown (i.e., the optimal phase for a 5 km vertical wavelength was searched in each case). In Figure 6 we show the SABER  profile. We used the same band-pass filter between 2 and 10 km as with the lidar data. Notice a clear nearly 5 km vertical periodicity. The sounding was located about 334 km northeast from Río Gallegos, so the 50 km minimum distance requirement was fulfilled. Equation 4 was used to calculate the missing component k * . The total horizontal wavelength that was found is 155 km and the deviation from the East direction was 2.2 • anticlockwise. The SABER horizontal displacement for measurements from 30 to 40 km height is 0.18 • northward and 0.17 • eastward, which is 23.5 km with an angle 58.7 • anticlockwise from the east direction. If this displacement is projected on the horizontal wave vector direction, it is equal to 13 km, which is small compared to the 155 km total horizontal wavelength. Therefore, in this case the satellite profile can be considered vertical. In order to quantify the possible distortion of the vertical wavelength observed in the slanted SABER profile, we use the formula derived by de la . According to the elevation angle of the sounding (23.1 • ) and the orientation angle of the wave vector found when it is projected on the vertical plane defined by the retrieval (86.7 • ), the possible error is about 12% (both angles are defined with respect to the ground). The GW intrinsic period was also calculated with our solution and the aid of the horizontal velocities provided by the ERA Interim reanalysis and was about 2.2 hr (19.6 m/s is the average speed parallel to the horizontal wave vector between the mountain tops around 600 mb, where presumably the waves are generated, and the maximum altitude of our study, at approximately 3 mb). The aspect ratio (the division of vertical and horizontal scales) and intrinsic period of this wave belong to the hydrostatic nonrotating regime, close to the border with the nonhydrostatic spectral sector (Gill, 1982).
Notice that changing the choice of the initial pair (one profile from the lidars and one from SABER) leads from Equations 3 and 4 to a different 2 × 2 linear equation system to be solved but both are equivalent, which means that they lead to the same solution. We now evaluate the impact of uncertainties in the specified parameters of any of both equation sets, which are represented by a 2 × 2 M matrix (four distances) and a column vector b (two phase differences). The equation set is then represented by and to evaluate if in our procedure any small changes in the known parameters can produce large changes in the solution s, we must calculate the condition number K(M). Then (e.g., Alexander & de la Torre, 2010) where s , M , and b refer to the relative error of s, M, and b, respectively. K = 1.47 for any of both equivalent equation sets in our example (1 is the optimal value in any case). Distances can be given with high precision, so their uncertainty can be estimated to be around 1% ( M = 0.01). Phase differences are extracted from the comparison of the same mode in the profiles at the two different places at the same time, and we could evaluate them to be around 10% ( b = 0.1). Then, variations in the solutions would be around 15%. If the three profiles would become nearly collinear, then the equation set would tend to an ill condition, and the propagation of the precision errors to the solution would have dramatic effects through the increase of K(M). In general, it is possible to calculate the horizontal phase speed c h and the horizontal and vertical components of phase velocity c x,y,z as seen from the ground for nonstationary waves The GW-associated specific potential energy can be also calculated over at least one wavelength as whereas the specific horizontal MF components may be obtained by an expression valid in the midfrequency range (e.g., Ern et al., 2004) where g is gravity, refers to the wavelengths, N can be obtained from any lidar profile, and finally T o and T are the GW temperature amplitude and the background atmospheric temperature (which can be obtained from either lidar), respectively. We obtain E p = 44.6 J/kg, F x = −1.45 J/kg, and F y = −0.06 J/kg. The average density in the studied height interval is 1.19 × 10 −2 kg/m 3 , and both components of MF then, respectively, become −0.018 and −0.001 Pa.

Conclusions
We have shown in a zone with high GW activity through one illustrative example that it is possible to reconstruct the 3-D structure of a dominant wave observed by two close simultaneous lidar soundings and an additional vertical temperature profile. The described method may help to dodge a present difficulty as is the determination of the three signed components of the wave vector. This is an essential quantity to determine directional MF, which is strongly related to atmospheric model parameterizations of GW drag. The last element is an Achilles heel that affects the simulation of zonal mean wind and temperature structure at middle and high latitudes. In particular, orographic waves like those generated in the hot spot here studied are currently considered to make significant contributions to the vertical transport of GW MF (e.g., McLandress et al., 2012).
Diverse tests verified if all instruments are likely observing the same GW and no important side effects contaminate the analysis. Only one out of 17 cases provided concrete results. This should be attributed to the fact that in several cases the diverse instruments were observing different phenomena or GW were too weak or affected by other effects or underwent diffusion, absorption, unstable, or nonlinear behavior. It should be noted that in some previous studies it was assumed that much longer latitude, longitude, and time intervals contained the same GW with no additional verification of coherency. McDonald (2012) evaluated as a function of horizontal separation and time difference the percentage of GPS radio occultation paired profiles that may be expected to contain the same GW. For example, his estimations were that approximately 30% of the pairs at 50-60 • S in the SH separated by less than 250 km and by less than 15 min were seeing the same GW. It is clear that an improved version of the Rio Gallegos lidar and eventually the 24 hr operation of both instruments would lead to a substantial increase in the number of cases meeting the required conditions.

10.1029/2020EA001074
The assumptions of the method should be recalled to constrain its validity: one dominant GW is observed by both lidars and the additional sounding, GW and background are adequately separated by the digital filter, and there are no further aliasing effects in the horizontal plane than those considered.

Appendix A
In Figure A1 we may see the spectral power for both lidars as a function of and s for the case that failed to meet the final test in order to be considered a possible standing GW. In Figure A2 we see the corresponding temperature perturbation against time and height for both lidars. The time frame of the coincident data for both sources is from 7 June 2018 01:04 (UTC) to 08:30 on the same day. It may be seen that in Rio Grande an upper stationary front of maxima is clearly defined for only approximately the last 1/3 of the total time.   Figure 5, it should be considered that the different aspect ratio is due to the fact that both cases do not have the same time extension.