Tropical and Subtropical South American Intraseasonal Variability: A Normal-Mode Approach

Abstract

Instead of using the traditional space-time Fourier analysis of filtered specific atmospheric fields, a normal-mode decomposition method was used to analyze South American intraseasonal variability (ISV). Intraseasonal variability was examined separately in the 30–90-day band, 20–30-day band, and 10–20-day band. The most characteristic structure in the intraseasonal time-scale, in the three bands, was the dipole-like convection between the South Atlantic Convergence Zone (SACZ) and the central-east South America (CESA) region. In the 30–90-day band, the convective and circulation patterns were modulated by the large-scale Madden–Julian oscillation (MJO). In the 20–30-day and 10–20-day bands, the convection structures were primarily controlled by extratropical Rossby wave trains. The normal-mode decomposition of reanalysis data based on 30–90-day, 20–30-day, and 10–20-day ISV showed that the tropospheric circulation and CESA–SACZ convective structure observed over South America were dominated by rotational modes (i.e., Rossby waves, mixed Rossby-gravity waves). A considerable portion of the 30–90-day ISV was also associated with the inertio-gravity (IGW) modes (e.g., Kelvin waves), mainly prevailing during the austral rainy season. The proposed decomposition methodology demonstrated that a realistic circulation can be reproduced, giving a powerful tool for diagnosing and studying the dynamics of waves and the interactions between them in terms of their ability to provide causal accounts of the features seen in observations.

1. Introduction

Intraseasonal variability (ISV) in convection over South America is associated with large-scale tropical, subtropical, and extratropical atmospheric disturbances and is manifested in different bands of the intraseasonal time scale: 10–20, 20–30, and 30–90 days, with different associated dynamics. For instance, Liebmann et al. [1] (see their Figure 2) and Paegle et al. [2] showed that there are distinct spectral peaks in these bands in the most prominent ISV over South America, in the South Atlantic convergence zone (SACZ) and central-east South America (CESA). Silverio and Grimm [3] also separated these three bands for southern Africa convection and connected them with corresponding bands of intraseasonal variability over South America, confirming findings by Grimm and Reason [4] on teleconnections between the two continents in the intraseasonal time scales. There was a peak around 12 days associated with quasi-biweekly oscillation [5], another peak around 22 days [2], and an MJO-related peak in the 30–90 day band (e.g., [2,6]). Liebmann et al. [1] found ISV peaks near 50, 27, and 11 days over the SACZ.

One of the most distinctive features characterizing the South American wet season (October–April) is the SACZ. It varies on many time scales, displaying prominent ISV in all its frequency bands. The SACZ variability is often part of a dipole-like configuration between CESA (where the strongest center of variability is found) and southeastern South America. It is stronger in the wet season (austral summer), in which it is the first ISV mode [6,7,8,9,10,11,12,13,14]. This mode (as well as SACZ) is primarily modulated by transient disturbances [1,2,6,8,13,15,16]. Rossby wave trains, which can be forced by tropical/subtropical convective activity, such as the Madden–Julian oscillation (MJO), induce ISV over South America [6,11,17,18,19]. This interaction between tropics and extratropics is frequently linked to the development of the Pacific–South America (PSA) teleconnection pattern (e.g., [20]). In submonthly time scales, Liebmann et al. [1] found two preferred paths of Rossby wave train patterns in the Southern Hemisphere: one affecting the SACZ and another influencing the southwestern Amazon. The southern Amazon pattern resembled the “cold surges” phenomenon discussed in detail by Garreaud and Wallace [21], Garreaud [22], Lupo et al. [23], among others. During the dry season (June to August), the ISV convective features were different from the wet season, especially in the 30–90 day band, when the strongest variability center moved southward. In the 10–30-day ISV over South America, a dipole-like structure (CESA–SACZ) was still visible, but with a much stronger signal over the CESA region during the dry season [11].

As described above, a large portion of intraseasonal rainfall variability is a result of the complex interaction in the 30–90-day ISV itself and shorter. In this line, Gelbrecht et al. [19] used a phase synchronization technique to corroborate that the CESA-SACZ dipole-like precipitation structure is caused mainly by extratropical Rossby waves [e.g., [7,17]]. However, some limiting factors of their approach include the irregular/intermittent character of phenomena often misrepresented by linear techniques such as empirical orthogonal functions (EOFs), as well as the lack of detailed attribution of the types and wavenumbers of the modes associated with the CESA–SACZ convection pattern. Meticulous treatment of the ISV over South America can be carried out using more intrinsically nonlinear approaches like self-organizing maps [SOM, [24]] rather than a traditional linear technique such as EOFs. Moreover, the wide variety of spatial and temporal variability in the intraseasonal time scale can also be addressed using normal-mode functions (NMF), which are orthogonal eigenfunctions of the linearized primitive equations on a sphere [25,26]. Indeed, recent works used NMFs to characterize physical properties representative of the MJO [27,28], and other tropical atmospheric disturbances [29,30,31]. An earlier study Baer [32] suggested a two-dimensional index ($index=s+n$) as a measure of horizontal scale, as in Kasahara [33] [see his Figure 5], where s and n are the zonal wavenumbers and meridional indices, respectively. Thus, in this study, we project 3D atmospheric fields onto normal modes of the global primitive equations, based on the Kasahara and Puri [25] approach, to determine the modes that more closely describe the observations.

Particular research issues relevant to intraseasonal oscillations over South America, which were not completely explored in previous studies, are addressed in the present study. We further explore the different mechanisms associated with the ISV over South America, which are associated with the aforementioned distinct spectral peaks, and present a multivariate three-dimensional analysis of the intraseasonal circulation based on normal mode expansion. Thus, to assess the physical mechanisms associated with the 30–90-day ISV, 20–30-day ISV, and the 10–20-day ISV, we computed a decomposition of these frequencies bands in terms of NMF [25] by performing linear regressions between the time series of the dominant mode at each season and time-scale and normal-mode amplitudes.

This paper is organized as follows. Section 2 presents a brief description of the data and methodologies. In Section 3.1 and Section 3.2, we described the observed dynamical mechanisms associated with South American intraseasonal variability. Section 3.3.1, Section 3.3.2 and Section 3.3.3 analyze the normal-mode components related to the South American 30–90-day ISV, 20–30-day ISV, and 10–20-day ISV, respectively. Finally, the main results are summarized and discussed in Section 4.

2. Materials and Methods

2.1. Satellite and Reanalysis ERA5

Satellite-observed OLR data were used as a proxy for convection. The OLR data were obtained from the National Oceanic and Atmospheric Administration—NOAA [34]. To assess the dynamical features and the modal decomposition of the ISV, we made use of daily data from the fifth reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF) [ERA5; [35]]. Both datasets cover the time period starting in 1980 and ending in 2016. The dataset has a horizontal spectral resolution of N80 (approximately ${1.125}^{\circ }$× ${1.125}^{\circ }$) (This denotes a lower resolution configuration of the ERA5 (32 km) TL639 grid) and 137 vertical levels ranging from 1012.04 up to 0.01 hPa. The variables used in this study were geopotential height (z), temperature (T), horizontal winds (u, v), specific humidity (q), and surface pressure ($sp$).

2.2. Filtering and Empirical Orthogonal Function (EOF) Technique

Daily anomalies in the convection and dynamical fields were calculated at each grid point by subtracting the first three harmonics (i.e., the annual cycle and two subsequent harmonics) for the entire 37-year time series, in order to remove the seasonal cycle. The intraseasonal filtered anomalies were obtained by applying fast Fourier transform (FFT), considering a frequency domain of 30–90 days, 20–30 days, and 10–20 days. Filtered OLR within the South American domain (${40}^{\circ }$ S–$5^{\circ }$ N and ${70}^{\circ }$ W–${30}^{\circ }$ W) was then submitted to a covariance matrix EOF analysis that retained the local variance in the EOF fluctuations. As in Kiladis et al. [36], EOFs were computed considering the entire record (from 1980 to 2016) but centered on each day of the calendar year using a sliding window. Then, 121-day, 91-day, and 61-day window lengths were considered for the 30–90-day band, 20–30-day band, and 20–10-day band, respectively. This approach considered the complex convective propagation over the region and better characterized the seasonal variation in the ISV. The time series of the standardized first principal component (PC1), for each frequency, were used to explore the main features associated with the respective time scales. The PC1 time series were also used in the normal-mode analysis.

2.3. Linear Regression

This study’s convective and circulation structures were based on linear regression. We regressed the standardized PC1 of the 30–90-day ISV, 20–30-day ISV, and 10–20-day ISV against dynamical and convective fields (OLR, velocity potential, streamfunction, and winds at 200-hPa). Day 0 of the composition was defined as the day with the lowest PC1. The statistical significance of these results was assessed based on a two-tailed Student’s t-test. This method takes into account the correlation coefficients and an effective number of independent samples (degrees of freedom) based on the decorrelation time-scale, as in Livezey and Chen [37] [more details in [38,39]]. For the sake of simplicity, we only considered two sub-seasons for the analysis: October–April, the period of active South American monsoon season, and May–September the period of the non-monsoon season.

2.4. Global Normal-Mode Function (NMF) Expansion

Given that the reanalysis data are provided for the entire globe, it was desirable to associate them with normal modes of the equations on the sphere. The linearized system of the atmospheric primitive equation in sigma coordinates in the vertical direction is given by

$$\begin{array}{c}\hfill \frac{\partial u^{\prime }}{\partial t}-2\Omega v^{\prime }sin\left(\varphi \right)=-\frac{g}{acos\left(\varphi \right)}\frac{\partial h^{\prime }}{\partial \lambda },\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\partial v^{\prime }}{\partial t}+2\Omega u^{\prime }sin\left(\varphi \right)=-\frac{g}{a}\frac{\partial h^{\prime }}{\partial \varphi },\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left[\frac{\partial }{\partial \sigma }\left(\frac{g\sigma }{R{\Gamma }_0}\frac{\partial h^{\prime }}{\partial \sigma }\right)\right]-\nabla ·{\mathbf{V}}^{\prime }=0,\end{array}$$

(3)

where $\varphi$ is the latitude, $\lambda$ is the longitude, ${\mathbf{V}}^{\prime }=(u^{\prime },v^{\prime })$ is the velocity field given by its zonal and meridional components,

$$\begin{array}{c}\hfill \nabla ·{\mathrm{V}}^{\prime }=\frac{1}{acos\varphi }\left(\frac{\partial u^{\prime }}{\partial \lambda }+\frac{\partial }{\partial \varphi }\left(v^{\prime }cos\varphi \right)\right).\end{array}$$

(4)

$\Omega$ is the Earth’s rotation rate, and a is its radius. $h^{\prime }=P/g$ represents the modified geopotential height, with P being the pressure field, g the acceleration of gravity, R and the gas constant for dry air, $\kappa =R/C_p=2/7$, where $C_p$ is the specific heat at constant pressure. ${\Gamma }_0=\kappa T_0/\sigma -d{T}_0/d\sigma$ is the static stability parameter, where $T_0=T_0\left(\sigma \right)$ is the globally horizontally averaged temperature. The boundary conditions are no-penetration conditions at the top and at the bottom ($\sigma ={\sigma }_T$ and $\sigma =1$).

The solutions of this coupled system were obtained by performing a separation of variables into a horizontal and a vertical structure:

$$\left[\begin{array}{c}{u}^{\prime }(\lambda ,\varphi ,\sigma ,t)\\ v^{\prime }(\lambda ,\varphi ,\sigma ,t)\\ h^{\prime }(\lambda ,\varphi ,\sigma ,t)\end{array}\right]=G\left(\sigma \right)\left[\begin{array}{c}u(\lambda ,\varphi ,t)\\ v(\lambda ,\varphi ,t)\\ h(\lambda ,\varphi ,t)\end{array}\right],$$

(5)

where the vertical structure function is given by $G\left(\sigma \right)$ and is expanded in terms of a basis of orthonormal basis functions:

$$G\left(\sigma \right)=\sum _{m=1}^M{c}_m{G}_m\left(\sigma \right),$$

(6)

where $G_m\left(\sigma \right)$ are the eigenfunctions of the vertical structure eigenproblem, such that

$${\int }_{{\sigma }_T}^1{G}_m\left(\sigma \right)G_n\left(\sigma \right)d\sigma ={\delta }_{mn},$$

(7)

where ${\delta }_{mn}$ is the Kronecker delta, which is equal to 1 if $m=n$, and equal to 0 otherwise. The coefficient $c_m$ is calculated by

$$\begin{array}{c}\hfill c_m={\int }_{{\sigma }_T}^{1}G\left(\sigma \right)G_m\left(\sigma \right)d\sigma ,\end{array}$$

(8)

${\sigma }_T$ is the model top in $\sigma$-coordinate. The horizontal structure-function is given by the product of an oscillatory term in time and a spatial structure, as follows:

$$\mathbf{U}(\lambda ,\varphi ,t)=\sum _{n=0}^N\sum _{k=0}^K{\mathbf{H}}_n^k(\lambda ,\varphi )e^{(-i{\omega }_n^{k}t)},$$

(9)

where n and k are the meridional and the zonal mode indices, respectively. The spatial structure ${\mathbf{H}}_n^k$ is described by Hough modes. Hough modes are categorized into Rossby–Haurwitz waves (RH), inertia–gravity waves (IG), mixed Rossby–gravity waves (MRG), and Kelvin waves (K). RH and MRG waves are rotational waves (ROT). The IGW modes are composed of K and IG waves. Given a vector field $\mathbf{X}$, on a discrete grid over the sphere, the projection of $\mathbf{X}$ onto the basis of normal mode function is obtained using the inner product from the vertical eigenvalue problem:

$$\mathbf{X}(\lambda ,\varphi ,{\sigma }_j)=\sum _{m=1}^M{\mathbf{X}}_m(\lambda ,\varphi )G_m\left({\sigma }_j\right),$$

(10)

providing a set of horizontal structures ${\mathbf{X}}_m$, for each vertical level $j=1,...,J$, and $G_m\left({\sigma }_j\right)$ is a discretized version of the vertical structure function $G_m\left(\sigma \right)$ via finite differences. ${\mathbf{X}}_m$ is then projected onto the basis of Hough functions to obtain the normal-mode coefficient associated with indices $(m,n,k)$:

$${\chi }_{mnk}={\int }_0^{2\pi }{\int }_{-\pi /2}^{\pi /2}{\mathbf{X}}_m·{\left[{\mathbf{H}}_n^k\right]}^{*}sin\left(\varphi \right)d\varphi d\lambda .$$

(11)

Based on this, the vector field $\mathbf{X}$ is expressed as a sum of components corresponding to each of the elements of the basis of the normal mode functions with their respective amplitude, and the index${ }^{*}$ represents the complex conjugate of the Hough mode:

$$\mathbf{X}(\lambda ,\varphi ,{\sigma }_j)=\sum _{m=1}^M\sum _{n=0}^N\sum _{k=0}^K{\chi }_{kmn}{G}_m\left({\sigma }_j\right){\mathbf{H}}_n^k(\lambda ,\varphi ).$$

(12)

In this study, we used the open-source software MODES [40], which performed these operations on the given ERA5 reanalysis. In other words, given a set of observed (reanalysis) horizontal winds and modified geopotential height fields evolving in time $\mathbf{W}={(u(\lambda ,\varphi ,\sigma ,t),v(\lambda ,\varphi ,\sigma ,t),h(\lambda ,\varphi ,\sigma ,t))}^T$. The modified geopotential height is a new variable that accounts for variations in the surface pressure around the globe, and it is defined as $P=\Phi +R{T}_{0}ln\left(p_s\right)$ [25], where $\Phi$ is the geopotential height, and R is the gas constant of the air. We used the discretized inner product (replacing the integrals by summations over the grid points) defined by the combination of Equations (8) and (11) to project the observed field onto the basis of normal mode functions. This provided a unique decomposition of the observed fields. Thus, the amplitude ${\chi }_{mnk}$ of the mode with zonal wavenumber k, meridional index n and vertical index m associated with $\mathbf{W}$ at time t is given by

$$\begin{array}{cc}\hfill {\chi }_{mnk}\left(t\right)& =\langle \mathbf{W}(\lambda ,\varphi ,\sigma ,t),{\mathbf{N}}_{mnk}(\lambda ,\varphi ,\sigma )\rangle \hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill & ={\int }_{{\sigma }_T}^1{\int }_0^{2\pi }{\int }_{-\pi /2}^{\pi /2}\mathbf{W}(\lambda ,\varphi ,\sigma ,t)·{\left[{\mathbf{H}}_n^k(\lambda ,\varphi )G_m\left(\sigma \right)\right]}^{*}sin\left(\varphi \right)d\varphi d\lambda d\sigma \hfill \end{array}$$

(13b)

where the normal mode function ${\mathbf{N}}_{mnk}$ is the product of the vertical structure function $G_m\left(\sigma \right)$ by the horizontal structure function ${\mathbf{H}}_n^k(\lambda ,\varphi )$.

One of the main advantages of this approach is to systematically attribute a certain type of atmospheric wave (i.e., Rossby, inertio-gravity, Kelvin, mixed Rossby-gravity) to observed fields. In this scenario, the role of nonlinear terms and the momentum and energy sources/sinks are included in the phase space-time evolution equation given by

$$\frac{d{\chi }_{mnk}}{dt}-i\omega {\chi }_{mnk}={\eta }_{mnk}+f_{mnk}$$

(14)

where ${\eta }_{mnk}+f_{mnk}$ represent the projection on the nonlinear terms and forcing terms in physical space on the normal mode ${\chi }_{mnk}$. Thus, the amplitude and phase of a particular normal mode (such as a Kelvin wave with the vertical structure given by $G_m\left(\sigma \right)$ with meridional mode n and zonal wavenumber k changes in time according to the impact of the combined effect of the nonlinearities and the physical forcing. ${\eta }_{mnk}$ represents the role in the interaction of all possible modes for ${\chi }_{mnk}$. Thus, in the absence of nonlinearities and forcing, the particular mode represented by the Kelvin wave should maintain its amplitude and the non-dispersive phase speed is the theoretical value that is approximate $\sqrt{g{h}_m}$, where $h_m$ is the eigenvalue of the vertical structure equation. Thus, Equation (14) shows that any deviation in the linear theoretical phase speed can only be attributed to the role of the nonlinearity and forcing. In a linearized state about a climatological zonal flow, the effect is included in both ${\eta }_{mnk}$ and $f_{mnk}$ if the basic state is not constant. The forcing term projection is required in order to have a stationary basic state.

In the usual interpretation of the Wheeler–Kiladis diagram [41], the reference Matsuno dispersion relation is provided in the background for a particular vertical mode with equivalent depth that more closely represents the influence of deep tropospheric diabatic convective heating ($\sqrt{g{h}_m}\cong 30$ ${\mathrm{ms}}^{-1}$). Therefore, when spectral energy is found along with the theoretical Kelvin regime, this means that free Kelvin waves contain a substantial amount of spectral power. However, the role of nonlinearities and forcing may distort the linear propagation speed (eventually inverting the direction) and cause substantial time change in the evolution of the Kelvin mode amplitude. Through the normal mode decomposition, we are able to detect spectral regions with significant distortion from the linear behavior caused either by nonlinearities (including the basic state role) and/or forcing.

However, there are some disadvantages to using this approach. For instance, the normal modes of the primitive equations are obtained through the linearization of a basic state at rest and ignoring physical processes such as radiative and diabatic processes and the presence of humidity [42,43]. Such processes might be important in the coupling between waves and convection [14,44,45,46]. In this scenario, the effect of the basic state on the atmospheric wave results from the nonlinear interaction between the waves and the basic states described as a superposition of normal mode functions. Furthermore, the choice of the basis of NMF is not unique but is a result of the chosen atmospheric model, and it has a clear physical interpretation, allowing us to associate particular observed atmospheric oscillations in terms of free-dry atmospheric waves.

2.5. Computation of the South American Intraseasonal Variability in Modal Space

In this work, we used indices that describe the tropical and extratropical precursors associated with high and low-frequency intraseasonal variability. The resulting precipitation pattern in South America was assessed in terms of NMF expansion following [27,30]. These authors introduced a methodology to decompose the widely-used all-season multivariate MJO index [47] into normal mode functions. They performed a linear regression of the RMM indexes $Y_i\left(t\right)$, $i=1,2$. The realtime multivariate MJO (RMM) index, is an index widely used to characterize MJO activity. It is based on a multivariate EOF analysis of zonal winds at 200 and 850 hPa and outgoing longwave radiation (OLR) in the equatorial Indo-Pacific region [47]. The index is then regressed against the normal mode function coefficients:

$${\mathcal{R}}_{kmn}^i=\frac{1}{N-1}\sum _{t=1}^N\frac{({\chi }_{kmn}\left(t\right)-{\overline{\chi }}_{kmn})(Y_i\left(t\right)-{\overline{Y}}_i)}{\mathrm{Var}\left(Y_i\left(t\right)\right)}$$

(15)

where ${\mathcal{R}}_{kmn}^i$ is the regression coefficient of the i-th index against the normal mode with mode indices $(k,m,n)$. In addition, the complex coefficient of ${\mathcal{R}}_{kmn}^i$ describes the projection of the Southern Hemisphere circulation associated with the 30–90-day ISV, 20–30-day ISV, and 10–20-day ISV. The relative importance of each normal mode to the given i-th index is obtained by its variance as follows:

$$V_{kmn}^i=g{D}_m{\mathcal{R}}_{kmn}^i{\left({\mathcal{R}}_{kmn}^i\right)}^{*}$$

(16)

where g is the acceleration of gravity, $D_m$ is the equivalent depth of the m-th vertical mode, and${ }^{*}$ indicates the complex conjugate, respectively (more details in [27]). In this work, the equivalent depth is the eigenvalue of the vertical structure problem, where each equivalent depth is associated with a different vertical profile. For example, the vertical mode $m=1$ has the largest equivalent depth $D_1=\mathrm{10,165.05}$ m and corresponds to the external mode (it does not change sign along the vertical coordinate). The vertical modes $m\ge 2$ are internal modes, and they change sign $m-1$ times along the vertical coordinate $\sigma$.

3. Results

3.1. EOF Results

Figure 1a shows that the intraseasonal time scale (10–90 days) over South America occurred in three main ranges. These modes occurred in pairs when the EOF analysis was computed for the entire 10–90-day band. The first and second modes (PC1 and PC2) peaked at about ∼48 days and together explained 26.8% of the total variance. The large amount of variance in PC1 and PC2 is typically associated with the large-scale MJO [48,49]. The second pair (PC3 and PC4) explained 14.5% of the total variance and peaks at about ∼15 days and ∼12 days. Similar spectral peaks were also documented in the OLR data for the submonthly time scale in Liebmann et al. [1]. Finally, the third pair (PC5 and PC6) peaked at about ∼24 days and together explained 9.3% of the total variance. Silverio and Grimm [3] found similar periods for convection over southern Africa, connected with that over South America, and consistent with the teleconnections documented in Grimm and Reason [4]. The spatial structure of the leading EOFs for the October–April and May–September seasons is presented in Figure S1 of the Supplementary Material. In the next section, the convection features and their associated circulation anomalies are investigated in detail for each band shown in Figure 1.

3.2. South America Intraseasonal Variability

3.2.1. South American 30–90-day ISV

Figure 2 shows the maps of OLR (shading), velocity potential ($\chi$; contours), and winds (vectors) at 200 hPa, obtained by regressing these fields onto the PC1 of the 30–90-day ISV. The lags considered in Figure 2 are based on the spectral peak (around 50 days) observed in Figure 1b.

Figure 2 depicts the evolution in tropical convection and implied large-scale circulation (upper-level velocity potential and winds) during a typical oscillation for the October–April (left column) and May–September seasons (right column). The large convection anomalies along the equator in the Maritime Continent (Figure 2a) propagated eastward to South America, resulting in a CESA–SACZ dipole-like configuration (Figure 2c). Despite the EOF calculations being made within the South American domain, the upper-level structure exhibited a zonal wavenumber-1 structure in the equatorial belt, as in other EOF-based analyses of the MJO. For instance, day 0 (Figure 2c) showed a strong positive center over the Maritime Continent and a negative center of action over South America.

On the other hand, Figure 2 (right column) showed a large-scale convection pattern and upper-level divergence for the May–September period. Some differences in eastward propagation phase speed for the circulation and OLR anomalies are evident in the regression maps. The OLR anomalies, initially over the Maritime Continent (day-25 in Figure 2a), propagated eastward creating conditions for convection over a broad area of South America (day-25; Figure 2c). The convective evolution, from day 25 on, clearly resembles the spatial structure widely documented using diverse MJO indices (e.g., [36,47]). There is an upper-level divergence center (negative velocity potential) over the Indian Ocean and a convergence center (positive velocity potential) over South America. The convective–dynamical evolution observed during October–April, even during the May–Sep months, clearly resembles the composite maps made using the OMI index in Mayta et al. [50]. In addition, during this period, the organization of the convective features was more strongly modulated by the Rossby wave trains (Figure A1), as well as being in agreement with the results presented in Vera et al. [11]. Because the PC1 of the 30–90-day ISV was calculated using a sliding window, our results demonstrated that this index properly represents the large-scale seasonal MJO impacts [51,52].

3.2.2. South American 20–30-Day ISV

Figure 3 shows the resulting lag regression horizontal composite based on the 20–30-day band. The October–April period, up until day-12 (Figure 3a), showed a wave train mainly centered over South America. This circulation pattern triggered enhanced convection over the CESA region. At day-4, (Figure 3b) the ridge located at about ${40}^{\circ }$ S propagated eastward due to a zonally oriented wave train propagating along that latitude. At the same time, enhanced convection over the SACZ appeared. Zonal westerlies were also triggered around the equator, resembling the mechanisms documented in Mayta et al. [39]. In addition, the 20–30-day band oscillations appeared to be linked to tropical convection resembling the extratropical and tropical teleconnections between South America and Africa documented in Silverio and Grimm [3], corroborating the teleconnections found in Grimm and Reason [4]. At day 0 (Figure 3c), enhanced convection over the SACZ region was observed as the wave train propagated equatorward. During the May–September period, a similar convective and circulation pattern was noticed. Compared to the October–April period, the wave train was more zonally oriented due to the strong winter jet stream.

3.2.3. South American 10–20-Day ISV

As in the previous section, here the 10–20-day ISV was regressed against dynamical and convective fields. Figure 4 shows the regressed values of OLR (shading), 200 hPa streamfunction (contours) for days 7, 4, and 0. The lags considered in Figure 4 are based on the spectral peak (around 15 days) observed in Figure 1b. During the Oct–Apr period, at day 7, as in Vera et al. [11], enhanced convection occurred over the CESA region (Figure 4a). At the same time, a well-developed series of upper-level alternating cyclones and anticyclones extending eastward and equatorward were observed. Then, 3 days later, suppressed convection signal started over Argentina, as the time of the wave trains propagated towards the South Atlantic Ocean (Figure 4b). At day 0 (Figure 4c), convection peaks over the SACZ region and the Rossby wave trains were propagating equatorward. Our results are consistent with previous works [1,2,53,54,55], who showed similar OLR and large-scale features associated with the submonthly variability over the SACZ region.

On the other hand, during the May–September season, the enhanced (day 7) and suppressed (day 0) convection covered a broad area of South America, with a northwest–southeast band extending to the adjacent ocean. Recently, Vera et al. [11] also showed similar strong convection over the CESA region, extending their signal towards the southern Amazon, as observed in Figure 4a. Similarly, the convective activity observed during the dry period was accompanied by highly statistically significant Rossby wave trains. The wave trains, unlike the Oct–Apr period, stretched eastward and equatorward from the Western Pacific and South Pacific Convergence Zone (SPCZ) with an arch-like structure. In addition, streamfunction anomalies were stronger in the dry season, centered at lower latitudes, and the OLR anomalies were much stronger over CESA than over SACZ. These wave trains also appeared to have a circumnavigating wave structure, with a clear wavenumber 6 in the Southern Hemisphere. This pattern resembled the spatial features associated with the PSA-like mode documented by Mo and Higgins [20]. A large portion of the extratropical wave trains in the upper-troposphere (Figure 4) were directly forced by the divergent outflow from regions of enhanced equatorial convection, such as the MJO convection [56,57]. The wave trains can be also forced by the subtropical compensatory convergence or even by subtropical convection, as shown in Grimm [6] regarding the extratropical teleconnection to South America that influences the observed circulation and convection patterns.

3.3. Normal-Mode Components of South American Intraseasonal Variability

The interaction between tropical convection and large-scale systems is characterized by energy conversion processes [58]. In this sense, the analysis of normal-mode decomposition of the intraseasonal variability constitutes a methodology for diagnosing the energy responsible for circulation. In this approach, the solution of the vertical structure associated with intraseasonal variability enables the analysis of the energetics for each of the vertical modes, separately, in external and internal modes. On the other hand, the energy distribution between the horizontal modes, the eigenvalues (normal modes) are classified into balanced modes (rotational, ROT) and unbalanced modes (inertio-gravity waves; IGW) (Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10, and Figure A2). The balanced modes include the Rossby wave (RH) and mixed Rossby-gravity (MRG) waves. The unbalanced modes include the Kelvin (K) and inertio-gravity (IG) waves.

3.3.1. Normal-Mode Decomposition: 30–90-day ISV

Figure 5 displays the contribution of each mode (zonal and vertical) to the total variance. At large scales (k = 1–5), most of the 30–90-day ISV variance was associated with ROT modes rather than IGW. Indeed, for k = 1, ROT alone accounted for 20% in the dry season and 50% in the wet season. Our results are in agreement with previous works [27,59,60], where the authors documented the same planetary modes for the MJO. We note that as the zonal wavenumber increases the spectral power of the inertio-gravity waves approaches that of the rotational waves, becoming equal at k∼30. A spectral cross over at k∼30 was previously reported in Žagar et al. [61]. Figure 5b,d reveal the leading vertical modes with a strong contribution of ROT modes for barotropic modes (Vertical modes with $1<m<5$ changed sign more than once along the vertical domain, but they only changed their sign at the stratosphere and kept a barotropic structure in the troposphere (see Figure 4a in Žagar et al. [40])) ($m<$ 5), while barotropic IGW modes were less prevalent. In this study, we determined a barotropic mode as every mode with an equivalent barotropic structure in the troposphere. For more internal modes ($m>$ 5), the contribution of IGW modes became more important, but they were still around one order of magnitude lower than the ROT modes. The distribution of energy spectra was similar throughout the year. The most noticeable difference was the larger contribution of Kelvin waves during the wet season.

Figure 6 shows the regression horizontal structure associated with the 30–90-day band for a pressure level close to ∼200–hPa. The projected circulation represents the contribution of the ROT modes (Figure 6a,d), IGW modes (Figure 6b,e), and the total fields (Figure 6c,f). The calculations were computed separately for the October–April and May–September seasons, respectively, and at lag 0 only (as in Figure 2c). Figure 6c,f suggest that a large portion of the 30–90-day ISV circulation, such as the mid-latitude wave-trains, was dominated by ROT modes (as in observation; see Figure A1). The ROT mode with k = 1 was the dominant mode associated with the tropical MJO [27], but in the extratropics, other low zonal numbers also made important contributions. These widely documented mid-latitude wave-trains (e.g., [6,11,19]; and references therein) present different aspects when comparing the dry and wet seasons, which can be explained in terms of the spectrum of their ROT variance in each season (Figure 5). Indeed, these wavetrains acquired a more clear pattern during the dry season, since in this season there was less energy for global scale wavenumbers $k=1-3$. This result was also expected, since the dominant 30–90-day ISV pattern from May–September was mainly influenced by extratropical disturbances rather than large-scale MJO eastward-propagation (Figure 2 and Figure A1). Figure 2b also depicts anomalous subsidence over South America, resembling the structure of the $k=1$ eastward propagating IGW mode (i.e., Kelvin wave). As observed in Figure 5c, for almost all k values, the contribution of the Kelvin waves was higher in the wet season than in the dry season. A relatively strong IGW signal over the Andes, as observed in Figure 6e, was a result of its interaction with the Southern Hemisphere winter upper-level westerlies that are stronger at this latitude. To better represent the MJO upper-level zonal wavenumber-1 (k = 1) structure in the equatorial belt, in Figure A2, we plotted the velocity potential instead of the streamfunction. The upper-level wind anomalies were mainly zonal (Figure A2b)with a wavenumber-1 structure comparable to those in previous studies [48,62,63]. The pattern was suggestive of an equatorial Kelvin wave signature that extended from South America, being barely equatorially trapped with a band of westerlies between ${10}^{\circ }$ N/S. Upper-level divergence over South America for the dry season, even in IGW modes, highlighted the presence of wave trains (Figure A2e).

3.3.2. Normal-Mode Decomposition: 20–30-day ISV

In this section, we analyze the normal mode contribution to the 20–30-day ISV over South America, which showed a spectral peak around 24 days (Figure 1b). As observed for the 30–90-day band, the variance was also dominated by rotational modes in both seasons (Figure 7). During the dry season, the largest contribution came from a zonal wavenumber $k=4$ for the ROT modes. The IGW contribution was of the order of 3–5% less than the ROT modes. Their contribution became relatively more relevant as k increased. For instance, for $k>20$, the contributions to the total variance from the IGW and ROT modes were almost equal. Kelvin waves made an important contribution (∼3%) to zonal wavenumber $k=1$, and its contribution decayed rapidly as k increased. The largest contribution to the MRG waves came from zonal wavenumbers k = 6–7. Similar variance spectra were observed during the wet season for the 20–30-day band (Figure 7c,d). The most important difference was observed in ROT modes, where the maximum contribution came for zonal wavenumber k = 1–2 instead.

The dominant terms for the vertical structure in both seasons, as expected, came from ROT modes with a barotropic structure in the troposphere ($m\le 5$). The contribution of baroclinic modes, on the other hand, was most relevant for m = 7–12.

The 20–30-day ISV regression patterns are shown in Figure 8 for a level close to 200 hPa. Figure 8c,f suggest that the majority of the average 20–30-day band circulation in both seasons was ROT. The intensity of the patterns was clearly stronger during May–September compared with October–April. The reconstructed circulation, as in observation (Figure 3), revealed a strong teleconnection between South America and southern Africa.

3.3.3. Normal Mode Decomposition: 10–20-day ISV

Following Equation (16), Figure 9 shows the contribution of the various modes to the 10–20-day ISV. The energy spectra were similar to the 20–30-day band, with the distribution of the variance dominated by ROT modes. The variance distribution on the zonal mode index k shows that the 10–20-day band was strongly dominated by ROT modes for large-scale modes (k = 1–7). On the other hand, as expected, the contribution of equatorially confined modes such as Kelvin and MRG modes was less relevant compared to their contribution to the 30–90-day ISV (Figure 9a,c).

The vertical distribution of the variance shows that the 10–20-day band was more associated with modes with barotropic structure in the troposphere (m = 1–5; Figure 9b,d). Considering that the 10–20-day band represents higher latitude dynamics, lower-order modes with the barotropic mode were dominant [25,64]. In addition, the contribution of modes with baroclinic structure was observed in modes with large m (m = 6–15), with peaks at m = 8–9. Peaks at m = 8–9 were quite evident for both seasons, and stronger during the May–September period, as also documented in Silva Dias and Bonatti [64]. The most significant difference was the large contribution of MRG modes for $m<$ 3 during the dry season (Figure 9b). Indeed, this could be explained by the fact that this mode has an asymmetric wind structure with respect to the equator and can have different responses owing to solar forcing depending on the time of the year [58].

The decomposition of the regressed circulation fields (upper-level streamfunction and winds) onto IGW and ROT components associated with the 10–20-day band is presented in Figure 10. According to Figure 10c,f, the average circulation features are given by ROT components, which was also expected from Figure 9. In other words, we could reconstruct the basic features of the previously observed structures in Figure 4 by using just ROT modes. The same predominance of ROT modes was found for both seasons. Compared with the circulation patterns of the 20–30-day ISV fields, during the dry season, it is noticeable that the IGW modes made a more important contribution.

4. Conclusions and Discussion

In this study, we presented an alternative approach to analyzing tropical and subtropical South American intraseasonal variability (ISV), based on normal-mode decomposition. This methodology, based on the 3-D multivariate decomposition in terms of the normal-mode functions, provides the relative contributions of the balanced (e.g., Rossby, mixed Rossby–gravity waves) and unbalanced (e.g., Kelvin, inertio-gravity waves) modes building the South American intraseasonal circulation. This information can be useful for studying the dynamics of waves and the interactions between them in terms of their ability to provide causal accounts of the observed features at different intraseasonal time scales. Furthermore, this method deals with global data, and thus, provides information on teleconnections in response to a determined forcing (e.g., [30]).

The leading patterns in the intraseasonal time scale were first obtained through EOF analysis, as in Kiladis et al. [36]. EOFs were computed onto the region of maximum intraseasonal variance. The period considered for the analysis was from 1980 to 2016, but centered on each day of the calendar year, using a sliding window approach to take into account the seasonal migration of the intraseasonal convection [36,52]. The results revealed three dominant patterns: (i) 30–90-day band, (ii) 20–30-day band, and (iii) 10–20-day band (Figure 1). Spectral analysis also confirmed significant peaks at ∼48 days, ∼24 days, and ∼15 days, respectively (Figure 1b). The most distinct feature, in the ISV for all bands, was the presence of a dipole-like CESA–SACZ structure, as documented in many previous studies [7,8,11,19]. However, the mechanisms responsible for the convective features were different in each corresponding band:

(i) In the 30–90-day ISV, the convective-circulation features were primarily caused by the large-scale eastward-propagating Madden–Julian oscillation (Figure 2). This modulation was more evident during the October–April period. During the May–September season, the resulting convection pattern (stronger over the CESA region) was mainly controlled by extratropical Rossby wave disturbances (Figure A1). During this period, a negligible modulation from the large-scale MJO was detected (Figure 2).

(ii) In the 20–30-day ISV, the wave trains and the corresponding convection structure over South America were likely more related to a direct forcing by the upper-level divergent outflow from equatorial convection (Figure 3). In addition, as documented by Liebmann et al. [65] and more recently by Silverio and Grimm [3], in this intraseasonal band there was a strong tropical and extratropical teleconnection between South America and southern Africa (day-4 in Figure 3).

(iii) In the 10–20-day ISV, our result showed that the CESA–SACZ convective structure was also modulated by the extratropical Rossby wave trains, as documented in many previous works [2,6,11]. While the 20–30-day ISV peaked at the submonthly time scale, the wave trains in the 10–20-day ISV were more synoptic-scale related. These wave trains also appeared to have a circumnavigating wave structure with a clear wavenumber 6 in the Southern Hemisphere (Figure 4). In contrast to the 20–30-day ISV, the extratropical wave trains were apparently not directly forced by the upper-level divergent outflow from enhanced equatorial convection, but by the subtropical compensatory convergence or even by subtropical convection, such as convection in the South Pacific convergence zone [6,17].

The relative importance of the balanced and unbalanced modes in the South American intraseasonal dynamical circulation patterns was also assessed in the present study. The idea behind this method is to seek the “flavors” responsible for observed circulation. Thus, by applying regression between the complex expansion coefficients of the NMF of the reanalysis data and time series of the 30–90-day ISV, 20–30-day ISV, and 10–20-day ISV, our results showed that rotational modes (e.g., Rossby waves, mixed Rossby–gravity waves) governed the tropospheric circulation in the three bands (Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10 and Figure A2). Other modes such as the inertio-gravity waves and Kelvin waves contributed considerably to the observed circulation in the 30–90-day band, mainly during the October–April period (Figure 5, Figure 6b and Figure A2b). Regarding the zonal and vertical structure, the larger contribution (for three bands) came from modes with the lowest zonal wavenumbers and the tropospheric equivalent barotropic modes (As depicted in Figure S2 of the Supplementary Material, vertical modes 1, 2, 3, 4, and 5 displayed barotropic characteristics in the troposphere. This implies that the vertical structure function of these modes remained consistent without changing sign throughout the troposphere. To clarify this distinction, we refered to vertical modes 1 through 5 as barotropic equivalent modes. Additional explanation and Figure S2 have been included in the Supplementary Material to provide further insights) (Figure 5, Figure 7 and Figure 9).

Subseasonal to intraseasonal variability over South America entails a complex and nonlinear interaction. In this study, we focused on the average modal picture of South American intraseasonal variability, while an alternative would be the NMF representation of selected events, such as the extreme drought event that occurred in the 2013/2014 austral summer [66]. The normal mode approach, therefore, serves as an alternative method for assessing intraseasonal variability in the region. The suggested decomposition methodology, distinguishing between low- and high-frequency intraseasonal components, can offer valuable insights into the dynamics of intraseasonal variability in South America. This approach presents an alternative tool for diagnosing model-related issues by comparing the normal mode decomposition of reanalysis data with model predictions of precipitation.