Investigation of the long-term variation of gravity waves over South America using empirical orthogonal function analysis

The spatial and temporal variability of gravity waves (GWs) potential energy (Ep) over South America (SA) was examined by analyzing temperature profiles obtained through the utilization of Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) from January 2002 to December 2021. We used the empirical orthogonal function (EOF) analysis to decompose GWs parameters and to analyze the GW variations over SA. We considered the first three eigenmodes (EOF1, EOF2, and EOF3) and their principal components (PC1, PC2, and PC3) of the EOF decomposition, which accounts for ∼ 80–90% of the total GWs variation over SA. Further, we analyzed the coupled variation of Ep and zonal mean wind ( U ) to verify their inter-dependencies using the singular value decomposition (SVD). The spatial variation showed that different localized mechanisms generate GWs at different sectors of the continent. The EOF1 of Ep comprised more than 50%, the EOF2 ∼ 20–25%, and the third ∼ 10–15% of the total GWs variation. The positive variation of GWs energy in the EOF1 is localized in the tropical region from the lower stratosphere to the lower mesosphere and southward below 1.5° S in the upper mesosphere. The spectral analysis of GWs energy showed biannual, annual, semiannual, and 11-year variations at different eigenvectors. Relative Ep (REp) showed an asymmetric hemispheric response to solar flux over South America. The REp response to QBO showed a modulating effect below 70 km and a positive response above 70 km. There is a good positive correlation between the temporal component of EOF2 of Ep and the quasi-biennial oscillation (QBO) at 30 mb and 50 mb in the PC2 temporal variation.


Graphical Abstract 1 Introduction
Atmospheric gravity waves (GWs) possess many parameters, including wave amplitudes, potential energy (Ep), momentum flux, vertical and horizontal wavelengths, frequencies and velocities, etc.An important parameter is the Ep, as the energy transport by breaking these waves at higher altitudes drives part of the mesosphere and lower thermosphere (MLT) as well as ionospheric mechanisms (Trinh et al. 2018).In recent decades, observations of GWs have contributed to a better understanding of dynamic saturation mechanisms, vertical propagation, and temporal and geographic variations (Fritts and Alexander 2003).In addition, computer modeling has revealed the characteristics and sizes of potential GW sources, the allowable spectral patterns in each layer of the atmosphere, energy transfer, and wave-to-wave interaction with the mean wind (Fritts and Alexander 2003).
The importance of GWs arises from their influence on atmospheric circulation and the thermal condition of the middle atmosphere (Vincent 1994).These are key reasons for improving the understanding of how global and regional activities result in the generation, propagation, and dissipation of GWs.A lot of research has been done over the years on the features of GWs using sensors on the ground (Dewan et al. 1998;Medeiros 2003;Smith et al. 2005;Suzuki et al. 2009;Li et al. 2011;Medeiros et al. 2018;Llamedo et al. 2019), satellite observations, and limb sounding methods (Dewan and Picard 1998;Tsuda et al. 2000Tsuda et al. , 2011;;Schmidt et al. 2005Schmidt et al. , 2010Schmidt et al. , 2016;;De la Torre et al. 2006, 2014;Alexander et al. 2008aAlexander et al. ,b, 2010;;Faber et al. 2013;Tsuda 2014;Ern et al. 2018Ern et al. , 2016)).Recent advancements in our comprehension of GWs have been significantly enhanced through a combination of in situ, ground-based, and space-based observations and theoretical and numerical investigations.The GWs significantly impact the overall circulation patterns of the middle atmosphere and lower thermosphere at a regional and global scale (Trinh et al. 2018).
The investigation of gravity waves has been conducted by Krebsbach and Preusse (2007) as well as Ern et al. (2011) utilizing the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument data.Krebsbach and Preusse (2007) conducted a spectral analysis on the time series of weekly zonal root mean square GW amplitudes using temperature data from SABER/TIMED from January 29, 2002, to January 2006.Their analysis focused on identifying systematic interannual and annual GW activity.In a separate study, Ern et al. (2011) investigated the absolute values of GW momentum flux in the stratosphere and the entire mesosphere using global temperature measurements.In this study, an alternative approach was employed to examine the SABER data spanning a period of 8 years.The objective is to statistically discern the primary wave sources, assess their spatial and temporal fluctuations, and investigate the fundamental attributes of the gravity waves originating from these sources.The breaking GWs have a significant role in the quasi-biennial oscillation (QBO) observed in the tropical lower stratosphere, as indicated by studies conducted by Ern and Preusse (2009) and Kawatani et al. (2010).Based on the temperature profile data, Liu et al. (2019) analyzed the trends in global activities of gravity waves from 2002 to 2015, as well as the relations of gravity wave activities to solar activity, the QBO, and the El Nino Southern Oscillation (ENSO).According to the SABER data, the GW activities in the altitude range from 20 to 100 km were highlighted.The regional long-term trends of GWs at altitudes between 20 and 100 km, especially over continents like South America (SA), with dramatic and various sources of GW yet to be thoroughly examined.
The application of empirical orthogonal functions (EOFs) in the field of meteorology was initially performed during the latter part of the 1940s.The methodology that involves the decomposition of a space-time field into spatial patterns and corresponding time indices has significantly contributed to advancing our understanding of the atmosphere.Nevertheless, because of the diverse range of atmospheric phenomena present, such as stationary and propagating elements, EOFs are limited in their ability to offer a comprehensive depiction.For instance, eigenfunctions are commonly challenging to interpret due to their geometric characteristics, overarching nature, and orthogonality in both spatial and temporal domains.Several modifications, such as the implementation of rotated EOFs, have been added to acquire more localized features.According to Hannachi et al. (2007), the EOF analysis is a good technique for analyzing atmospheric parameters.This technique has also been used in various studies before, e.g., (Matsuo and Forbes 2010;Flynn et al. 2018).
The primary objective of this study is to examine the long-term variations of GW distribution in the stratosphere and mesosphere (20-100 km) over SA.The SA continent encompasses different localized meteorological mechanisms that generate GW in different sectors, which have been relatively understudied in the context of long-term investigations.The remainder of this paper is organized as follows: the methodology section describes the technique, which includes the datasets and methods for extracting GWs and GWs parameter decomposition.The result section examines GW's spatial and temporal variations.The discussion section will explain the physical characteristics of the results before the conclusion section.

Methodology
We used continuous TIMED/SABER temperature profile data from January 2002 to December 2021.The GWs characteristics were extracted from satellite-observed temperature profiles (T) using scale separation.The main concept is to separate the small-scale perturbations ( T ′ ) from the background temperature ( T ).The T ′ profile or the post-processed (e.g., band-pass filtering or spectral decomposition and reconstruction) GWs properties can be determined (Preusse et al. 2002;Ern et al. 2004Ern et al. , 2011Ern et al. , 2018;;Alexander et al. 2008a;Liu et al. 2014).To calculate the GWs parameters, first, we found the mean of all the raw temperature profiles within every 10 o × 10 o latitude by longitude cell in all the selected pairs to get an effective background temperature (T) (Faber et al. 2013;Ayorinde et al. 2023).We used the maximum overlap discrete wavelet transform (MODWT) multiresolution analysis (MRA) (Percival 2008) to decompose the mean temperature profiles in a cell.MRA is the process of decomposing a signal into constituent parts that, when combined, provide the original signal perfectly.The decomposition of the signal is crucial for its usefulness in data analysis.The wavelet MRA separates the signal components using fixed functions known as wavelets or MODWT.Utilizing scaled waveform, MODWT may effectively detect local non-periodic patterns and signal singularities and define signal structures by measuring signal fluctuations while concurrently assessing the signal's temporal and scaling features.To extract patterns from the data, wavelet filtering is thus more effective than the conventional linear transport frequency filters (Percival and Walden 2000;Scafetta and West 2005).To extract the GWs from the temperature profiles measured by satellites, MRA wavelet decomposition was used to decompose the raw temperature profile ( T raw ) to extract the background temperature (T) as shown below: where T res is the small-scale temperature fluctuations caused by gravity wave activities, respectively.The primary goal is to distinguish between large-scale waves (such as planetary waves and tides) and small-scale fluctuations (T res ), as well as background temperature (Liu et al. 2022).Further, to remove noise and to derive the vertical wavelength ( v ), we applied the continuous wavelet transform (CWT) (Torrence and Compo 1998) to each T res , given as Then we applied the inverse CWT to each T ′(h, v ) and restricted the v to a range of 5-30 km to derive the temperature fluctuation ( T ′ ), remove noise, and obtain an altitude dependent v .Restricting the vertical wavelength will also enable us to remove waves that are not GWs (Kogure et al. 2017;Schoon and Zulicke 2018).We estimated the T amp using the T ′ and its Hilbert trans- form ( H(T ′) ) following (Kogure et al. 2017;Schoon and Zulicke 2018;Liu et al. 2022) and is given as: The Ep is therefore estimated using Eq.(4) below, (1) (3) where g is the acceleration due to gravity, N is the Brunt- Väisälä frequency, and T , and T ′ are background temper- ature and the temperature fluctuations caused by gravity wave activities, respectively.T amp in Eq. ( 4) is the same as Eq. ( 3), which is the temperature amplitude.The Ep calculation is based on the accurate extraction of T ′ in Eq. ( 2), and N is given as follows: where C p is the specific heat capacity of air at constant pressure, and h is the altitude.The gravity wave Ep is calculated using Eq. ( 4).The Ep depends only on the raw temperature profile, which can be separated into the background temperature (T), and the temperature fluctuation ( T ′ ).In calculating Ep, the T ′ is the variable needing the most careful attention.Figure 1a shows the raw temperature profiles taken within the altitude range of interest (from the stratosphere to lower thermosphere (20-100 km)).The raw temperature profiles are interpolated to 200 m intervals.The fluctuation profiles (Fig. 1b) are obtained by subtracting the background temperature profile from the raw temperature profiles.The temperature amplitude profile and the corresponding vertical wavelength (Fig. 1c, d) are obtained from the CWT.It is worth noting that the largest vertical wavelengths are observed between 40 and 70 km.The absolute Ep profile (Fig. 1e) is obtained using Eq. ( 4) for every temperature profile. (5) 3 The GW variation evaluation We analyzed the variations of GW parameters (Ep and Tˆ); first, we employed the empirical orthogonal function (EOF) analysis, which provides the spatial and temporal variability at different modes of variability to the total variance.Furthermore, it provides the expansion coefficients of the contributions of each mode in a single geophysical variable.To compute the EOF, we detrended the GW datasets along the third dimension (along the monthly data dimension) of the GW matrix using the learning least squares detrending following the method of Greene et al. (2019) to remove linear trend along axis from data, and later, we computed the EOF eigenmodes, principal components (PCs) of the EOFs, and their expansion coefficients.We employed the Eigen decomposition theorem to derive the EOF eigenmodes.The EOF analysis aims to identify a limited set of independent variables (predictors and factors) that capture the maximum amount of original information without any redundancy.EOF analysis is a useful tool for objectively examining the structure of variability in a data set, such as gravity wave Ep.It may also be used to evaluate correlations among different variables.The EOF analysis reveals the spatial patterns of the primary factors that contribute to the changes in gravity waves throughout time.Our primary aim was to construct a 2D matrix of the varying GW parameter (Ep) to identify coherent quantities extracted that are relevant to the behavior of the GW parameters across the latitudes and longitudes.The 2D matrix of the leading eigenmodes that can be used to reconstruct the GW datasets was represented as a 2D matrix of latitude and longitude, as described in Eq. ( 6).The reconstructed 2D matrix of Ep produces the spatial variations (eigenvalues) and also the temporal variations (eigenvectors).The EOF can represent the GW variability efficiently.This efficiency is attributed to the orthonormality of the eigenvectors and the magnitude ordering of the eigenvalues, which indicates the significance of each expansion term.
where Π denotes the Ep or T , EOF k signifies an orthogo- nal function independent of time, while Amp k represents the time-dependent PC amplitude of each EOF eigenmode.The index k ranged from 1 to the maximum number of eigenmodes contained in the GW function, with the eigenmodes ordered by decreasing variance.The variable d = 1,2,3,…, corresponded to the number of months between January 2002 and December 2021, and x (lati- tude and longitude) denoted the spatial positions associated with the observations.Remarkably, the first three EOF eigenmodes accounted for over 90% of the total variance and exhibited statistical significance.In contrast, the residual EOF eigenmodes contributed only a negligible difference, which was statistically insignificant.The distinct significance of the first three eigenmodes became evident due to their well-separated eigenvalues.As a result, our study focused primarily on these crucial three eigenmodes.
Also, to facilitate a deeper analysis of the GW parameters, we employ singular value decomposition (SVD) analysis to uncover the interconnected patterns of covariability among the GW parameters.The SVD and the EOF share a commonality in their focus on decomposing a covariance matrix.The SVD technique is very similar to the EOF method, but it excels in detecting the coupled modes of variability between geophysical variables and their temporal variations, thereby offering valuable insights into the percentage of covariance explained by each paired mode (Bretherton et al. 1992;Bjornsson and Venegas 1997).The SVD decomposes a matrix into three separate matrices, which can be useful for dimensionality reduction, noise reduction, and extracting important features from data.The SVD technique involves decomposing a "cross-covariance" matrix obtained from two variables.Mathematically, for a given matrix A, SVD represents it as the product of three matrices: ( 6) where A is the original matrix, U is an orthogonal matrix containing the left singular vectors of A. Σ is a diagonal matrix containing the singular values of A, and V T (trans- pose of V) is an orthogonal matrix containing the right singular vectors of A. In our study, the matrix U is the Ep, and the matrix V is the mean zonal wind U with dimen- sion of months of the year (in our case, 12 months for 20 years) by the multiplication of the size of latitude and longitude.Therefore, matrix A represents the matrix of the product of Ep as the left-eigen matrix and U as the right-eigen matrix.The results of these techniques (EOF and SVD) are presented in the next section.

Method for extracting long-term changes of GW
The relative Ep (REp) zonal-mean for each height level at the altitude range of 20 to 100 km were derived by averaging all the REp profiles within a single month and binned into 10° ×5° longitude-latitude grids, with a latitude overlap of 2. 5 o .We defined the REp as REp = Ep/Ep where Ep is the time-mean of Ep at each altitude.The linear trend of REp zonal-mean time series and the responses to solar activity (F 10.7 ), QBO wind, and ENSO are calculated for each latitude band and height level using the multivariate linear regression (MLR) approach (Wolter and Timlin 2011;Li et al. 2013).The MLR can be used to determine the link between one dependent variable (e.g., REp) and two or more independent variables (e. g F 10.7 , QBO wind, and ENSO).The MLR method's equation is as follows: The monthly zonal-mean value of GW parameters at month (j) and year (i) is represented by Π(t i,j ), and the quantity µ represents a constant k h value.The param- eter β 0 , which depicts the change in GW parameters over time, represents the monthly zonal-mean Ep linear trend from 2002 to 2021.The parameters β 1 , β 2 , β 3 , and β 4 , show the relationship between the time series of GW parameters and the time series of the four indices, depicting the monthly GW parameter's zonal-mean response to QBO at 30 hPa and 50 hPa, ENSO, and solar activity.The residual of the regression model can be utilized to estimate the standard deviation and p-value of each coefficient.This estimation can be achieved using the variance-covariance matrix and the student t test (Kutner et al. 2004;Mitchell et al. 2015).From 2002 to 2021, Fig. 2 shows the (8) Withi = 2002Withi = , 2003Withi = , ..., 2021;;and, j = 1,2, ..., 12 reference time series of QBO, ENSO, and solar activity.In Fig. 2a, the reference time series for 30 hPa is the temporal fluctuation of QBO zonal-mean zonal wind across the equator (Gavrilov et al. 2002).The QBO data are from radiosondes (obtained from the Meteorological Service Singapore Upper Air Observatory (1.34 o N, 103.89 o E), NASA (National Aeronautics and Space Administration)/GMAO (Global Modeling and Assimilation Office) assimilated data, and NASA satellites through the NASA earthdata site.The bi-monthly Multivariate ENSO index (MEI) values in Fig. 2b show the temporal variations of the ENSO phases (Wolter and Timlin 2011).In Fig. 2c, the solar activity is represented by the monthly mean values of radio emissions from the sun at a wavelength of 10.7 cm (F 10.7 ) (Geller et al. 2016), and F 10.7 is provided in solar flux units (sfu, where 1 sfu = 10 −22 Wm −2 Hz −1 ).We note that during the examined period, solar activity was at its peak from around 2002 to 2005 and 2012 to 2016, and lowest around 2007 to 2010 and 2017 to 2021.

Results
In this section, we present the results based on the analysis described in the methodology section.The 20-year mean GW parameters, the GWs variation through the decomposition method, and the MLR analysis investigated are presented.As seen in Eq. ( 4), the Ep is directly proportional to the temperature amplitude ( T ).There- fore, from Fig. 3, the representation of these two parameters (Ep and T ) are similar.Hence, for the rest of this study, our focus will be on the Ep only.In this section, we refer ± 20 o as the tropical region, 2 0 o S-40 o S as the sub- tropical region, and 40°S-55°S as the extratropical region.

The mean GW E p
Figure 3 shows the total mean of Ep in longitude and latitude over SA for the period of 20 years at four different altitude ranges (20-30 km and 30-40 km in the stratosphere, and 60-70 km and 80-90 km in the mesosphere).Also, 60-70 km and 80-90 km were chosen in the mesosphere to study the stratosphere-mesosphere coupling mechanism.GWs are seen to dominate the tropical region (± 20 o ) at 20-30 km as shown in the Ep values, extending to about 30°S.
At 30-40 km, Ep peaks over three major parts of SA: at 20° N, 20° S around the Andes Mountains, and over the Argentine and Patagonian mountains.The result at ∼20°N is suggested to be due to the northern hemisphere summer intertropical convergence zone (ITCZ), which generated gravity waves through convective systems.The two later results are suggested to be the GWs generated from the mountain waves and the influence of the jet streams from the polar region (Vadas and Becker 2019).We suggest that this phenomenon may not be obvious at 20-30 km because most of the low frequency generated around the equator have not encountered the mean zonal wind (Alexander et al. 2008c).From 30 km to somewhat around 70 km, the GW energy experiences damping because of the likelihood that most of these waves around the equator flow in the same direction as the eastward zonal wind, thereby reaching a critical level.
At 60-70 km and 80-90 km, the GW energy rise is seen to be centered around the southern Andes Mountain and over the Patagonian mountains, while the tropical GWs energy over SA is seen to be conserved with and increase around the equator at 80-90 km.We observe the lowest GW energy over the land area of central Brazil and Bolivia and closer to central America in the northern hemisphere.From Fig. 3, the GW energy difference over SA is ∼3.5-17J kg −1 with the larger energy seen at 20-30 km (∼17 J kg −1 ) and the lowest at 30-40 km (∼3.5 J kg −1 ).This result (the same for T ˆ (∼2 dB), figure not shown) clearly suggests the interaction of GWs at 30-40 km with the wind that could allow the wave to reach the critical level, thereby allowing the wave energy to be conserved.

Spatial and temporal variability of GW E p
The spatial variability of Ep is shown in Fig. 4 and its temporal variability is shown Fig. 5 at 20-30 km, 30-40 km, 60-70 km, and 80-90 km altitude ranges.As previously mentioned, we emphasized that the T EOF results will not be presented further since Fig. 3 and the EOF analysis (not shown here) has shown that the Ep and T are directly proportional and therefore, the results are apparently similar.The EOFs shown in Fig. 4a-d, illustrate spatially static patterns.Our spatial data was normalized between − 1 and 1.This guarantees uniform scaling of all features.If the original data has a broad range of values, this will help preserve the information contained in it.The values are dimensionless and normalized so that a value of 0 indicates the absence of variability in that region.The signs of the EOFs are initially arbitrary until they are assigned to their respective signed amplitudes in time (shown in Fig. 5).The initial three eigenmodes in the data set explain 92%, 91%, 95%, and 90% of the overall variability at 20-30 km, 30-40 km, 60-70 km, and 80-90 km, respectively.Due to the small amount of the explained variability associated with the rest of the EOFs, further discussion of these components will be omitted.We considered the yaw cycle of the TIMED/SABER satellite that cause a poleward discontinuity of the measurements at 52.5 o S-52.5 o N every 3 months.Therefore, we replaced the negligible space in the southern SA with the latitudinal average to minimize any statistical error caused by the yaw cycle.
At 20-30 km, the EOF1 accounts for about 52% of the total variation.Most of the positive variation is centered around ± 10 o apparently from the western Pacific decreasing towards the east.This phenomenon has also been reported by Ayorinde et al. (2023).This means that over 50% of GWs over SA at 20-30 km occur in this region.Deep convective activities characterize this area.It is also seen at ∼15 o S-40 o S (over the Andes Mountain range also often referred to as the subtropical region) that about 50% of the GWs attained zero variation.This means the mountain waves over the Andes at this height are vertically propagating waves which can extend up to  In the mesosphere at 60-70 km, EOF1 accounts for about 74% (Fig. 4) of the total variation of GWs over SA.Positive variation of GWs in the tropical region (about 20°N-25°S) showed the continuous propagation of GWs up to ∼70 km.EOF2 and EOF3 accounts for about 12% and 10%, respectively, of the total variation of GWs over SA.The positive variation of GWs is observed to be increasing from ∼60 km at ∼15 o S-40 o S (subtropical range) to ∼80 km extending to the mid-latitude with an increasing percentage (the result for the 70-80 km altitude range is not shown here).We also note the 15% variation in the EOFs at 40-50 km (result not shown) showed zero/negative variation.This finding is consistent with the research conducted by Liu et al. (2019), which revealed a pronounced GW peak above the Andes mountains during the Southern Hemisphere winter at ∼40 o S. The GW activity extends up to an altitude of ∼55 km.Additionally, it was discovered that orographic gravity waves exhibit a breaking phenomenon occurring at altitudes surpassing the maximum height of the stratospheric jet.At ∼55-65 km, the breaking and deposition of mountain waves (MWs) result in the production of body forces that give rise to secondary GWs on a larger scale (Bossert et al. 2017).We also observed resurgence of positive variation in EOF3 in the mid-latitude.At 80-90 km, the EOF1, EOF2, and EOF3 have ∼51%, ∼25%, and ∼14% of the total GWs variation respectively.Positive variation in EOF1 could be seen at ∼15 o S-55 o S in contrast to observation from lower altitudes.
The temporal variation of Ep is shown in Fig. 5 at 20-30 km, 30-40 km, 60-70 km, and 80-90 km altitude ranges.The PC1, PC2, and PC3 are the amplitudes of the EOF1, EOF2, and EOF3, respectively.We applied a linear polynomial fit on the estimated amplitude to quantitatively extract the patterns, understand the characteristics, and to identify important features of temporal variation.Also, we applied spectral analysis on the temporal variation to identify different variations in each eigenvector.At 2030 km, the spectral analysis in PC1 showed minor biannual and annual variations, PC2 had biannual, annual, semiannual, and 11-year variations, and PC3 had annual and semiannual variations.The temporal variation showed the positive peak in summer (December, January, and February) and the negative peak in winter (June, July, and August) with sometimes irregular amplitudes.For example, we observed the lowest positive peak in summer 2003, 2007, 2009, 2012, 2014, and 2021 in PC1 in Fig. 5d.These dates coincide with the sudden stratospheric warming (SSW) events ranging from minor to extreme SSW events when probably the QBO phase is westward (Siddiqui et al. 2018;Li et al. 2023).
The highest positive peak was found in January 2013 when there was an extreme and a trail cooling (TC) SSW event with an eastward QBO phase (Siddiqui et al. 2018;Li et al. 2023).The TC expresses the trailing upper stratosphere cooling anomaly strength with the trail phase duration more than 21 days (Li et al. 2021(Li et al. , 2023)).This result will be investigated in our further studies.The temporal variation at 30-40 km showed that PC1 had annual variation, PC2 had biannual, annual, and semiannual variations, and PC3 had annual and semiannual variations.Also, at 60-70 km, PC1 and PC2 showed annual variations, and PC3 had annual and semiannual variations.This depicts a stability in the GWs propagation in the upper stratosphere and lower mesosphere.Lastly, at 80-90 km, PC1 and PC2 showed semiannual and annual variations with a minor 11-year variation at PC2, and PC3 also showed annual, semiannual, and 11-year variations.The 11-year variation seen in PC2 and PC3 suggests the extent of solar cycle influence in the upper mesosphere and in the lower thermosphere.
The interesting result is the amplitudes in the EOF2 mainly at 20-30 km and 30-40 km.The amplitudes show a form resembling a QBO with a dramatic fluctuation at the peaks.We note that the percentage of EOF2 at 20-30 km and 30-40 km is ∼25% and ∼35%, respectively, which could be considered as substantial percentage.This result will be investigated further in Fig. 8. Above 40 km the amplitudes showed interannual oscillations with much smaller values up to ∼70 km.This could be as a result of the GWs interaction with the wind.PC3 showed terannual with interannual variation with much lower amplitudes.This could result from the GWs interactions with the local body forces (Vadas and Fritts 2001;Vadas et al. 2023).We noted that the correlation between the amplitudes and their respective fit showed correlation coefficients of ∼70-93%.The PC1 at 80-90 km has the lowest In Fig. 7, we employed the SVD technique to compare the coupled variations of the Ep and zonal wind parameter ( U ) using the MERRA-2 (Modern-Era Retrospec- tive analysis for Research and Applications, Version 2) zonal wind data.The SVD and EOF share similarities in their approach of decomposing a covariance matrix.
In the context of EOF analysis, the covariance matrix is constructed using a singular spatio-temporal field.The selection of the SVD approach is motivated by its ability to effectively capture patterns characterized by high covariance between two variables.Previous research has demonstrated its capability to adequately represent atmospheric and oceanic processes (Wilks 2015).The utilization of this approach is considered to be highly effective in examining the prevailing patterns of interaction, Fig. 7 The spatial and the temporal variation of the first and second mode of the normalized SVD analysis between Ep and U over SA over the 20-year period.The third column are the amplitudes of the normalized first and second squared covariance (SC1 and SC2).a is the altitude range of 30-40 km for SC2, b is the altitude range of 30-40 km for SC1, c is the altitude range of 20-30 km for SC2, and d is the altitude range of 20-30 km for SC1.The third column showed the amplitudes as it facilitates a better comprehension of the interconnections among various sets of variables (Frankignoul et al. 2011).
The coupled variation of the U and Ep using the SVD technique over SA is shown in Fig. 7 over the 20 year period.In Fig. 7, U was taken from the MERRA-2 Fig. 8 The temporal variation of the EOF2 at 30-40 km (a) and 20-30 km (b) over the 20-year period, and the QBO measurement at 30 mb (solid line in magenta) and at 50 mb (dotted line in red).c is the correlation coefficients (in %) of the Ep and QBO at 30 mb and 50 mb reanalysis data.The first mode of the coupled variation (SFC ∼97% variation) at 20-30 km and 30-40 km from Fig. 7d and c, respectively, showed that the spatial variation was in good agreement with positive variation in the tropical region at 20-30 km and in the mid-latitude at 30-40 km, a negative variation in the mid-latitude at 20-30 km and in the tropical region at 30-40 km.The temporal variation corresponds to a semiannual variation with the Ep with slight interannual variation at the peaks.However, as pointed out previously, the temporal variation in Fig. 7 is predominantly annual.There is only one peak each year.There is a good correlation coefficient between the Ep's and U 's (R = 67% at 20-30 km and R = 83% at 30-40 km) temporal variations.
The second mode of the coupled variation (SFC ∼97% variation) at 20-30 km and 30-40 km from Fig. 7c and Fig. 7a, respectively, showed a contrasting spatial variation with positive/negative Ep and U .The temporal vari- ation has a good correlation coefficient between the Ep and U (R = 66% at 20-30 km and R = 52% at 30-40 km).Interestingly, the temporal variation showed a quasibiennial and interannual variation like the PC2 result from Fig. 5c, d.The Ep's interannual variation could be because of the QBO on GWs.
To investigate further the quasi-biennial variation seen in PC2 in Fig. 5c and d and in SC2 in Fig. 7c and a, we took the PC2 from the Ep at 20-30 km and 30-40 km overlaid with the QBO measurement at 30 mb and at 50 mb shown in Fig. 8b and a at 20-30 km and 30-40 km, respectively.Also, Fig. 8c showed the correlation of Ep and QBO at 30 mb and 50 mb and QBO at 30 mb and 50 mb for all the altitude ranges.From Fig. 8a and  b, the second mode of the EOF (PC2) and QBO showed a correlation coefficient of ∼20% and ∼60% at 20-30 km and 40-50 km, respectively for 30 mb QBO.For 50 mb QBO, the correlation coefficient of ∼64% and ∼67% at 20-30 km and 30-40 km, respectively.The rest of the altitude ranges have negative or less significant correlations.This showed that the positive correlation of the Ep in PC2 with the QBO increases as the altitude increases in the stratosphere.We suggest that the second eigenvector of the EOF decomposition of a global gravity wave Ep or temperature variance needs to be evaluated for the prospect of the general relationship between the QBO index and the temperature in the stratosphere.
The correlation results showed that the second mode decomposition of Ep or U could be a very good alter- native to QBO indices.The QBO is characterized by zonally symmetric easterly and westerly wind regimes alternating regularly with periods varying from about 24-30 months (with an average of about 28 months) in the equator.QBO propagates vertically (downward) at a rate of 1 km/month.The oscillation is symmetric about the Equator with the average maximum amplitude of 20 ms −1 .Successive regimes first appear at the height of 10 mb and drop down to the height of 100 mb.The maximum amplitude is about 40 to 50 ms −1 conspicuously observed at altitude of 20 mb.The amplitude of the easterly phase is about twice as strong as that of the westerly phase (E = 30-35 ms −1 , W = 15-20 ms −1 ).The downward motion of the easterlies is usually more irregular than that of the westerlies.Transition between westerly and easterly wind is often delayed between 30 and 50 mb.Below 50 mb, QBO signal changes rapidly (i.e., rapid attenuation occurs below 23 km) (Holton 2004).

Ep trends and responses to ENSO, QBO, and solar activity
Figure 9a displays the deseasonalized trend of REp (% per decade), obtained using MLR fitting on the deseasonalized Ep data.Based on the data presented in Fig. 9a and f, it is evident that there exist both positive and negative trends in Ep, contingent upon the latitude and height variables.Negative variation in Ep (Fig. 9a) trends are found in the equatorial region extending to ∼50 km and in the higher latitudes in the lower stratosphere over SA.A positive trend could be seen in the northern SA stratosphere.These negative and positive trends could be attributed to strong planetary waves activities and SSW that cause the GWs modulations.During these events, the activity of GWs is significantly inhibited (Ern et al. 2016).The tropical and subtropical stratosphere is strongly affected by winds with different phases of QBO (Zhang et al. 2012).The above-mentioned apparent features exhibit a near symmetry across both hemispheres because of a sixmonth disparity in atmospheric conditions between the Southern and Northern Hemispheres (Ern et al. 2016;Liu et al. 2017).There are negative trends with the QBO at 50 mb at ∼50-75 km around the equator (Fig. 9e).This will be discussed further in the next section.
There is a substantial positive variation around ∼55-65 km near the latitude of ∼20 o S-40 o S which could be attributed to the positive trend in QBO at 30 mb and 50 mb.Ern et al. 2016 investigated the interconnections among GWs, SSW occurrences, and PJO events by examining various satellite data sources from 2001 to 2014.They demonstrated that, in regions of higher latitudes in the northern hemisphere, SSWs have the potential to modulate GW activity.The outcome is contingent upon specific variables, such as a polar vortex split and PJO occurrence.The potential influence of SSWs and PJO events may contribute to the observed positive trend.It could also be linked to the GWs breaking and momentum deposition.According to Liu et al. (2019), in the Southern Hemisphere at the approximate latitude of 40 • S, a notable concentration of GWs were observed above the Andes mountains.This concentration extends to an altitude of approximately ∼55 km.Based on their analysis of wind patterns and topographic data, these orographic GWs exhibit a breaking phenomenon occurring at altitudes surpassing the peak height of the stratospheric jet.At the altitudes ranging from ∼55-65 km, the dissipation of these GWs and the subsequent exchange of momentum result in the production of local forces that contribute to the emergence of more substantial secondary GWs on a wider spatial extent.
Figure 9b illustrates the relationship between the response of Ep and ENSO, as measured by the Multivariate ENSO Index (MEI).The analysis of Fig. 9b reveals a statistically significant positive/negative correlation between the response of Ep and ENSO within the latitude range of ∼20 o S-40 o S and at altitudes below ∼25 km.Gavrilov et al. (2004) showed that there exists a negative correlation between the intensity of GW seen by MF radar in Hawaii and the Southern Oscillation Index (SOI), which serves as a measure of Pacific ENSO activity.The SOI exhibits a negative correlation with El Nino events, which are characterized by warm conditions in the central Pacific.Conversely, a positive correlation is shown between the SOI and La Nina events, which are associated with chilly conditions in the central Pacific.Hence, the observed inverse relationship between the SOI and the intensity of GWs suggests that GWs exhibit greater intensity during El Nino episodes.In a similar vein, the observed correlation between the positive reaction of the REp and the MEI, as depicted in Fig. 6d, indicates that a robust El Nino event (characterized by a large positive MEI index) is associated with heightened intensity of the Madden-Julian Oscillation (MJO) and its associated convective activity.Therefore, the observed correlation between the reaction of the REp to the ENSO in the subtropical stratosphere aligns with the findings reported by Gavrilov et al. (2004).
The response of Ep to QBO at 30 mb (24 km) and at 50 mb (21 km) are shown, respectively, at Fig. 9c and d.We found a lesser or negative response of GW REp to QBO at 30 mb below 30 km and a positive response around 40 km at 0 o -20 o S. A high positive response at ∼20 o S-50 o S below ∼25 km followed by a negative response at around 30 km and 40 km.The QBO at 30 mb and 50 mb showed a relative positive response to REp at ∼20 o S-50 o S above 40 km.There is a high positive response of GW REp to QBO at 50 mb below 40 km at 0 o -20 o S and a high negative response at ∼20 o S-50 o S. In the equatorial region, there is a high negative response to the QBO at 50 mb above 40 to ∼70 km.Our result showed that the QBO can modulate the GWs energy at different regions of the SA atmosphere.Although the tropical stratosphere being the source of the QBO, its impact on GW activity can extend to higher latitudes in the mesosphere (Baldwin et al. 2001;Ern et al. 2008).It is observed that there is a general positive response of the REp to the QBO at approximately above 70 km over SA.The negative response observed in the stratosphere is in contrast to the opposite response observed in the mesosphere, since the phase of the QBO in the mesosphere is around 180 out of phase compared to the stratosphere (Liu et al. 2017).
Figure 9e illustrates the relationship between the response of Ep and solar flux (F 10.7cm ).We found two distinct responses of the REp to the solar flux: a general positive response in the southern hemisphere SA at ∼20-100 and a negative response at ∼30-90 km in the northern hemisphere of SA.This result is consistent with the result of Liu et al. (2017).Ern et al. (2011) demonstrated a negative association between GW momentum flux and solar flux, utilizing 9-year datasets from SABER at altitudes of 30-70 km.The observed decrease in REp in response to solar flux aligns with the intensity of GWs as detected by MF radar scans conducted in Saskatoon (52 o N, 107 o W) Gavrilov et al. (1995).Gavrilov et al. (2002) reported a positive association between GWs and solar activity.The relationship between GWs and solar activity appears to be influenced by various factors, such as the characteristics of the GW sources (such as nonuniform orography in the longitudinal direction) and the conditions of wave propagation across different longitudes (Gavrilov et al. 2002).

Discussion
The mean variations of GWs Ep from 20 years of SABER satellite data showed the morphology of GWs coupling between the stratosphere and mesosphere.GW energy is the amount of energy GWs propagates through the atmosphere, and the amplitude measured the maximum displacement of these waves as the propagates through the atmospheres.The relationships between Ep and T are the direct proportionality (Ep ∝ T 2 ) and the energy trans- ported increases as the amplitude increases.The concentration of GWs energy at ± 20 o at 20-30 km shows that GWs in the lower stratosphere are mainly generated from deep convective sources.Deep convection is commonly associated with this region (Liu and Vadas 2013;Llamedo et al. 2019;Ayorinde et al. 2023).Zhang et al. (2012) presented observational evidence to substantiate their claim that deep convection has a significant role in the reported activity of tropical gravity waves.Also, Ayorinde et al. (2023) found a good correlation between the tropical stratospheric GWs and the precipitable water vapor over SA consolidating the assertion that the tropical GWs are mostly generated from the deep convention.The effect of the tropical wind could be seen at 30-40 km by depleting the wave energy in the equatorial region.Kim et al. (2003) demonstrated that convective GWs above 20 km with phase speed (c) >> 0 ascend vertically until they break and reach critical level.The GWs generated convectively propagating in the same direction as the zonal wind reach their critical level in the middle stratosphere (30-40 km) but the GWs generated convectively propagating in the opposite direction as the zonal wind dissipated in the mesosphere (John and Kumar 2012;Liu et al. 2022).Other middle stratospheric characteristics (e.g., stratospheric radiative balance, ozone concentration, Brewer-Dobson circulation, wind patterns etc.) are also factors to be considered.
The results EOF1 (Fig. 4 first column) showed that different mechanisms are responsible for the generation and propagation of GWs dependently and independently over SA.Our result showed that regions with prominent excitation of GWs through convective system account for over 50% (On the average) of the total GWs variations.This means that GWs are mainly driven from convective sources, and they could propagate up to 70 km or 80 km before experiencing dissipation.This also depends on the season and other atmospheric conditions.For example, Zhang et al. (2012); Ern et al. (2016) posited that tropical deep convection serves as a significant contributor to the generation of gravity waves within the stratosphere.However, it has been discovered that the characteristics of gravity waves in the tropical/subtropical stratosphere are notably influenced by the presence of winds associated with various QBO phases (Li et al. 2023).
The mountains waves, jet mechanisms, coupled with convective systems mostly generated by frontal systems are likely to be responsible for the 25 to 35% of the total GWs variation (Fig. 4 second column) at 2 0 o S-40 o S. McLandress et al. (2000) and Jiang et al. (2002) analyzed data from the Microwave Limb Sounder (MLS) and deduced that the prominent changes in brightness recorded over the Andes Range can be attributed to the influence of mountain waves.Over this region also, Hierro et al. (2018) observed both orographic and convective activity above the Alps and the Andes mountains from radio occultation measurements, and in the Brazilian sector, Nyassor et al. (2022) observed a concentric GWs from an OH All-sky airglow imager at 29.44°S, 53.82°W attributed to a moving mesoscale system.The negative variation at 50-60 km agrees with result of Liu et al. (2019), indicating that GWs over the Andes could breakdown at this range and deposit momentum that generates local body forces, creating larger GW that later propagates to higher altitudes.Figure 4 (third column) accounted for 10 to 15% of the total variation which centered around the GWs generated by jet streams and the Patagonian mountains (Xu et al. 2018).
PC1 clearly indicates semiannual variation in the lower altitudes and terannual variations in the higher altitudes coupled with interannual variations of Ep.The semiannual variations are due to the tropical/subtropical convectively generated GWs in the summertime and the polar jet and mountain GWs generated in the wintertime.The terannual cycle could potentially be associated with the length of vigorous convection in subtropical latitudes.Krebsbach and Preusse (2007) observed a pronounced semiannual oscillation in the upper stratosphere/mesosphere at high latitudes.The annual components of the data exhibit two distinct peaks located at mid-latitudes in both the northern and southern hemispheres.These peaks are closely linked to the winter polar vortex.Additionally, there are peaks observed throughout the summer months in the subtropical regions.The indications for QBO are observed in the expansion of latitudes towards the mid-latitudes in the stratosphere of both hemispheres, as well as in the equatorial mesopause.The unprecedented correlation between the QBO and the Ep (evident in the T ) in Fig. 8 is also evident in the SC2 of the coupled variation of Ep and U in Fig. 7.
This occurrence of the biennial amplitudes in the stratosphere PC2 aligns with the descending westerly shear phase and can be ascribed to the QBO.The presence of a QBO signal in the GW activity has been detected through measurements of Global Positioning System Radio Occultation (GPS-RO) at latitudes around the equator within the lower stratosphere (Wu et al. 2006;De la Torre et al. 2006).Also, Liu et al. (2017) revealed that the impact of the QBO on the GW Ep is negative in the upper stratosphere of tropical regions with a positive response in the subtropical region.Furthermore, this negative response spreads to higher latitudes as altitude increases.The formation of the QBO is thought to be attributed to the deposition of momentum flux caused by radiative damping and critical level filtering by vertically propagating waves of different sizes (Antonita et al. 2008).Both gravity waves and equatorial planetary scale wave modes play significant roles in driving the QBO in the stratosphere (Ern et al. 2008).

Conclusion
The variabilities of GW activity in the stratosphere and mesosphere over South America were demonstrated based on latitude, longitude, and year.The estimation of GW parameters has been conducted using the temperature profiles acquired from the TIMED/SABER measurements spanning from January 2002 to December 2021.For the first time, we used the EOF and SVD to analyze 20 years of gravity waves E p over SA.In this study, we presented the first three eigenmodes (EOF1, EOF2, and EOF3) of the EOF and their respective principal components (PC1, PC2, and PC3).
GWs are seen to dominate the tropical region (± 20 o ) at 20-30 km as shown in the Ep.At 30-40 km, Ep is large in three sectors of SA: at ∼2 0 o N, ∼2 0 o S around the Andes Mountains, and over the Argentine and Patagonian mountains.At 60-70 km and 80-90 km, the GW energy rise is seen to be centered around the southern Andes Mountain and over the Patagonian mountains, while the tropical GWs energy over SA is seen to be conserved with and increases around the equator at 80-90 km.
Our study showed that the generation mechanisms of GWs accounted for more than 50% of the total variations of GWs over SA, over 25% are driven by large scale oscillations like QBO, and lesser than 20% are driven by other atmospheric phenomena.The temporal variation showed the positive peak in summer and the negative peak in winter with sometimes irregular amplitudes.The interesting result is the amplitudes in the EOF2 mainly at 20-30 km and 30-40 km.The amplitudes show a form resembling a QBO with a dramatic fluctuation at the peaks.The temporal variation is a semiannual variation with the Ep showing interannual variations.There is a good correlation between the Ep and U 's temporal variations.
Among the new findings, our result showed a very good agreement between the second mode of the EOF and SVD analysis with correlation of ∼20% at 20-30 km and ∼60% at 20-30 km for 30 mb QBO and correlation of ∼64% at 20-30 km and ∼67% at 20-30 km for 50 mb QBO.Tˆ and QBO showed a higher correlation of ∼34% and ∼79% for 30 mb QBO at 20-30 km.Also, Ep at 40-50 km showed a signature of 30 mb QBO with ∼42% correlation.While the rest of the altitude ranges have negative or less significant correlations.Also, there is a feasibility in the use of EOFs to understand more of the distribution and variation of GWs regionally and globally.The eigenmodes are linked to the generation, variation and mechanisms aiding or inhibiting the GWs propagation in the stratosphere and mesosphere.The SVD is a good technique to understand the inter dependencies of GW parameters.
Our analysis also includes examining the impact of solar activity, the QBO, and the ENSO on the GWs.The global response of REp to the solar flux is predominantly positive in the southern SA and negative in the northern hemisphere SA at 30-70 km.The response of GWs to the eastward phase of the QBO is characterized by both negative and positive responses below 70 km and positive above 70 km.The REp response to ENSO showed a positive response at 2 0 o S-40 o S at 20-25 km and a negative response above 25 km.

Fig. 1
Fig. 1 Example for data analysis using TIMED/SABER measurement to calculate the Ep and its parameters.The profiles are taken on 11 April 2019, a Temperature profiles at 20-100 km in altitude, b Temperature fluctuation from the temperature profile, c Derived temperature amplitude, d Vertical wavelength from the CWT analysis from the temperature fluctuation profile, and e The estimated Ep profile

Fig. 2
Fig. 2 Reference time series from years 2002 to 2021 were used for the regression analyses.a 30 hPa (24 km, blue) and 50 hPa (21 km, red) zonal winds over the equator to characterize the QBO (top panel), b Multivariate ENSO Index (MEI) to describe the ENSO signal (middle panel), and c Solar radio flux at 10.7 cm (F 10.7 ) to represent the solar flux (bottom panel)

Fig. 3
Fig. 3 The mean Ep over SA at four different altitudes during the period of 20-year period.a 60-70 km, b 80-90 km, c 20-30 km, and d 30-40 km.The Ep are measured in Jkg −1

Fig. 4
Fig. 4 The spatial variation of the eigenmodes of the normalized empirical orthogonal function (EOFs) analysis of the GW Ep over SA.The left column showed the first EOF eigenmode (EOF1), middle column showed the second EOF eigenmode (EOF2), and the right column showed the third EOF eigenmode (EOF3).a is the altitude range of 80-90 km, b is the altitude range of 60-70 km, c is the altitude range of 30-40 km, and d is the altitude range of 20-30 km.The color bars are dimensionless

Fig. 5
Fig. 5 The temporal variation of the eigenmodes of the normalized principal components (PCs) of the EOF analysis of Fig. 4 of the GW Ep over SA.The left column is PC1, middle column is PC2, and the right column is PC3. a is the altitude range of 80-90 km, b is the altitude range of 60-70 km, c is the altitude range of 30-40 km, and d is the altitude range of 20-30 km.σ represents the percentage of the covariance.The light blue lines are the linear polynomial fit of each of the temporal variation of the eigenmodes of the normalized PCs.σ represents the percentage of the covariance

Fig. 6
Fig. 6 Spectral analysis of the eigenvectors of the normalized principal components (PCs) of the EOF analysis of Fig. 5 of the GW Ep over SA.The left column is PC1, middle column is PC2, and the right column is PC3. a is the altitude range of 80-90 km, b is the altitude range of 60-70 km, c is the altitude range of 30-40 km, and d is the altitude range of 20-30 km.σ represents the percentage of the covariance

Fig. 9
Fig. 9 The deseasonalized linear trend term of relative Ep (a), along with the responses of relative Ep (b to e) to ENSO (MEI index), QBO at 30 mb, QBO at 50 mb, and solar flux (F 10.7 cm ), respectively examined in latitude-height sections.The black crossed lines are the regression coefficients where p-values are greater than 0.2.The regression coefficients with p values are less than or equal to 0.1 are denoted by the red lines and unmated.