The role of nonlinear self-interaction in the dynamics of planetary-scale atmospheric fluctuations

A central role in the general circulation of the atmosphere is played by planetary-scale inertial fluctuations with zonal wavenumber in the range k = 1–4. Geopotential variance in this range is markedly non-gaussian and a great fraction of it is non-propagating, in contrast with the normal distribution of amplitudes and the basically propagating character of fluctuations in the baroclinic range (3 < k < 15). While a wave dispersion relationship can be identified in the baroclinic range, no clear relationship between time and space scales emerges in the ultra-long regime (k < 5, period >10 days). We investigate the hypothesis that nonlinear self-interaction of planetary waves influences the mobility (and, therefore, the dispersion) of ultra-long planetary fluctuations. By means of a perturbation expansion of the barotropic vorticity equation we derive a minimal analytic description of the impact of self-nonlinearity on mobility and we show that this is responsible for a correction term to phase speed, with the prevalent effect of slowing down the propagation of waves. The intensity of nonlinear self-interaction is shown to increase with the complexity of the flow, depending on both its zonal and meridional modulations. Reanalysis data of geopotential height and zonal wind are analysed in order to test the effect of self-nonlinearity on observed planetary flows.


Introduction: the dispersion relation of atmospheric global scale fluctuations and the role of planetary waves self-nonlinearity
Intermittence and persistence are substantial components of atmospheric phenomenology: observations show the occasional appearance, in both the extra-tropical and tropical atmospheric circulation, of long-lasting features associated with specific and recurrent flow patterns-like blocking highs in middle latitudes and monsoons in the tropics-that are characterized, from a dynamical point of view, by enhanced stability with respect to ordinary flows. The idea that such persistent features may be more predictable found an explicit formulation in the 1979 work by Charney and DeVore [1] (hereafter referred to as CDV) in connection with blocking, and in Charney and Shukla [2] in connection with monsoons, raising great interest in view of the development of extended range forecast. Intensive research work following the seminal CDV paper highlighted a number of properties of the atmospheric low frequency variability (LFV) in the range 10-40 days: the non-normal, possibly bimodal [3,4], statistical behavior of mid-latitudinal geopotential fluctuations, as opposed to the unambiguously normal behavior of fluctuations in the short range (3-5 days) that, as well known, are essentially associated with the classical jet cycle (baroclinic instability-barotropic convergence), and the problematic behavior of average zonal wind (mean westerly momentum) which is subject to normally distributed fluctuations confined to a very restricted range [3,5].
Additional information concerning atmospheric LFV comes from space-time analysis of planetary scale fluctuations: since Blackmon it is known that systematic connections exist among space and time scales and the mobility of patterns of middle-latitude atmospheric circulation [6,7]. In fact, simple low pass filtering at 10 days allows one to separate of nontraveling global scale features from shorter-scale traveling fluctuations. All the above phenomenological knowledge finds a synthetic quantitative representation in Fourier spacetime Hayashi analysis [8,9] in the form of spectra revealing the existence in the ultra-long spectral range (zonal wavenumber k<5, period >10 days) of a quite robust non-traveling flow component [10]. Use of Hayashi spectra as global metrics in order to measure the ability of global models in representing the internal dynamical processes of atmospheric circulation shows that the dynamical nature of LFV still has to be fully understood in order to adequately model atmospheric circulation [11].
It is evident from the analysis of Hayashi spectra that a wave dispersion relationalthough complicated by baroclinic instability and other processes (in particular barotropic convergence)-holds at medium-small scales (k>4, corresponding at 50°of latitude to 8000 km), while no distinct dispersion relation can be identified in the ultra-long regime (k<4) where Rossby wave dynamics should be operating. It is well known that a number of relevant linear dynamical processes are operating in the ultra-long range which can account for the presence of a conspicuous fraction of non-propagating variance: flow modulation by mountains [12], heat sources [13], and, in particular, the presence around k=3 of a zero phase speed of the Rossby wave dispersion relation (the linear resonance of planetary waves). However, the effect of such processes does not seem to fully explain the observed phenomenology, in particular for what concerns the evolution of large amplitude waves. The point of view we propose in this work is that a key role in determining the observed phenomenology in the ultra-long spectral range is played by self-nonlinearity, i.e., the nonlinear interaction of planetary waves with themselves occuring when the amplitude of fluctuations is large. This process, local in wavenumber, can be seen as a possible alternative to the energy cascade process, i.e., the interaction of long waves with much shorter ones, initially proposed by Green [14] and further investigated by, e.g., Shutts [15] and Malguzzi [16], to explain the maintenance of blocking episodes. In order to explain the limited range of variations of zonal wind statistics and the intermittent behaviour of planetary waves, selfnonlinearity was already introduced in many theoretical studies as a possible mechanism capable of producing multiple equilibration of stationary Rossby waves in the presence of orographic relieves for a specified meridional profile of the zonal wind [17][18][19][20]. In this work, this approach is generalized to traveling planetary disturbances, confined in the extra-tropical atmosphere by the gradient of mean potential vorticity.
The paper is organized as follows. In section 2 we introduce free Rossby wave selfinteraction by means of a weakly nonlinear expansion in terms of wave amplitude of the barotropic vorticity equation. In sections 3 and 4 we look for the observational signature of planetary scale fluctuations self-interaction by analysing the ECMWF ERA-Interim reanalysis [21], considering the extended winter seasons October-April from 1979 to 2015. The selfinteraction coefficient defined by the theory is computed from the geopotential height of the northern hemisphere extra-tropics and, as suggested by theoretical treatment, a statistical correlation with the observed planetary wave phase speed is sought for. In section 5 conclusions are drawn.

Self-interaction of atmospheric planetary waves: a weakly nonlinear theory
Based on past analyses of weakly nonlinear interaction of Rossby waves through selfadvection [17][18][19][20], we assume that part of tropospheric LFV can be associated with selfinteraction of planetary fluctuations. This assumption is tested by computing the intensity of the nonlinear self-interaction of Rossby waves in the context of a barotropic framework. This kind of problem has been widely explored both for the atmosphere and the ocean: see, for example, the recent papers [22] and [23] and literature therein. Our objective is to derive from theory a quantitative formulation of the nonlinear self-interaction of planetary scale atmospheric fluctuations that is suitable for observational analysis. Mirroring the approach of several previous works, the adopted procedure is that of multiple-scale analysis. We do not restrict our study to the stationary case, thus considering any zonally propagating wave that is confined in the meridional domain (i.e., not crossing the equator). We derive an analytic relationship linking zonal phase speed to wave amplitude that can be possibly tested by looking for a statistical correlation in long time series of geopotential height from reanalysis data.
In order to model the mobility of planetary geopotential fluctuations and the effect of their self-interaction we need to effectively represent the advection of potential vorticity, which is the basic propagation mechanism of vertically coherent disturbances (as most of the planetary scale variance) in the troposphere. To this aim, we adopt the barotropic vorticity Its solution is sought in the form thus separating the symmetric component (corresponding to the zonal wind profile) from zonally asymmetric fluctuations. The expressions above are both dimensionless and the classic quasi-geostrophic scaling is assumed: the reference quantities are horizontal length, time scale and horizontal velocity taken equal to 10 6 m, 10 5 s, and 10 ms −1 , respectively. Spherical coordinates are employed, with λ and f indicating longitude and latitude. The Laplace operator is defined as a 2 cos cos cos 2 2 the jacobian the Coriolis parameter f = W f 2 sin , the rescaled Earth radius = a 6.37, and the dimensionless Earth rotation rate p W = / 2 0.864. The solution of the barotropic vorticity equation is not to be sought analytically, being preferable to construct a uniformly valid approximation of l f Y( ) t , , through multiple-scale analysis. As highlighted, for example, by Bender and Orszag [24] this is a particularly useful technique for the estimation of solutions of perturbation problems. The perturbation solution is obtained by expanding the streamfunction as a power series of ε in the form Consistently with classical Rossby scaling hypothesis, the zonal wind component corresponds to the zero order and the nonsymmetric part to higher order terms. The assumption that the latter component is small compared to the former directly entails that the order of magnitude of the perturbation parameter is given by the ratio of the amplitudes of the eddy and symmetric velocities, i.e.
In the adopted approach, the latitudinal wind profile f ( ) u is considered constant in time, clearly in contrast with its observed variability. Fixing the value of the time derivative of zonal wind as zero has the advantage of leading to a great simplification of the problem in discussion. In computing self-nonlinearity, the input zonal wind profile will be defined through a proper average in time. This choice is supported by the classical interpretation of the general circulation theory, revolving around the idea that turbulent fluxes have the effect of maintaining an 'equilibrium circulation', i.e., a time-mean circulation remarkably different from the true stationary solutions. Many dynamical theories in last decades have been founded on the hypothesis of a correspondence between time-mean and stationary configurations, although the naive application of averaging procedures might be misleading. These notions were discussed by [25,26]. The consistency of the classical theoretical framework was tested in [25] by analysing a model atmosphere on the basis of an 'equilibrium' circulation as approximated by the average zonal flow. A linear stability analysis of the modeled circulation proved it to be unstable. The conclusion was that several statistically independent processes contribute to the definition of a most probable state for circulation which carries no particular dynamical meaning.
In the past, however, LFV was also described as a consequence of the eddy-instability of a basic state constructed by perturbing the climatological mean flow with stochastic fluctuations, as discussed in [27]. Other authors such as [28] adopted in their work a forcing capable of keeping the climatological steady state in the absence of any disturbance and not acting on disturbances themselves. Likewise, we assimilate the input zonal wind to its time mean average, thus implicitly assuming the existence of a forcing allowing the zonal component of the solution to keep its steady state and not acting on the higher order asymmetric terms.
An approximate solution l f Y( ) t , , is found through multiple-scale analysis, which is applied to the present problem by introducing a new set of slow time variables {τ n } with n=1, 2, K such that t e = t. n n Even though the physical solution will be a function of time alone, in the employed formal procedure we will treat the new variables t and {ε n } as independent. We assume a perturbation expansion of the form u a cos and 3 t a n 2 . 7 2 The multiple-scale time derivative ¶ | t M is computed recurring to the chain rule for partial differentiation: so that, inserting the perturbation series (5) into (6) and collecting like powers of e, we obtain the following set of equations: , , 9 we also impose that the τ 1 derivative be zero, so that A=A(τ 2 ). Hence, the solution of (9c) assumes the form: where f ( ) p satisfies the following differential equation: The necessary and sufficient condition for the solvability of the third order equation (9d) is that the projection of the rhs on the solutions of the adjoint of (9b) be zero. Assuming the standard definition of scalar product in spherical coordinates with f cos as area element, the operator L is self-adjoint and the adjoint of equation (12) becomes 0 ,whose eigenfunctions are * = ¢ - Therefore, the condition to avoid secularity becomes: The above integral can be rewritten in the compact form with the coefficient χ (hereafter referred to as self-interaction coefficient) given by  [20], which was obtained in the case of stationary Rossby waves forced by orography. The solution of (16) can be written in the form: Since the simultaneous presence of the (arbitrary) coefficient ε and of the (arbitrary) amplitude A 0 in the perturbative solution becomes obviously redundant, we can set A 0 =1 without losing generality. Therefore, ε assumes the role of the dimensionless wave amplitude. Expression (18) is a correction term for the phase speed c. The solution for y 1 can be rewritten as: In the above expression of a second order correction to the phase speed is introduced which is directly proportional to the self-nonlinearity coefficient χ that, in turn, is an implicit function of the zonal wavenumber.

Numerical implementation
We estimate the strength of nonlinear self-interaction working on 6-hourly data from ECMWF ERA-Interim reanalysis [21]  The zonal average of the zonal wind field f ( ) u t , is computed and the index Δu(t 06 ) representing the variation of zonal wind between 12 and 00 UTC is introduced: with t 00 , t 06 and t 12 indicating each triplet of successive time steps at 00, 06 and 12 UTC. As our theoretical estimation of the strength of nonlinear self-coupling relies on the assumption of a zonal wind constant in time, we choose to compute χ at 06 UTC only from those time steps characterized by a small value of Δu(t 06 ). All 06 UTC time steps satisfying to this criterion are fed into the eigenvalue problem (12). A numerical solution of (12) is sought by reducing it to a matrix form and employing a finite difference scheme, which allows to compute the eigenvectors f { ( )} g and the associated eigenvalues { } c . We restrict our analysis to the first two eigenvectors, hereafter named modes A and B, as they are physically more relevant. Mode A is characterized by zero nodes, while mode B has zero or one node. All eigenvectors are taken with maximum amplitude equal to +1. Singular solutions are excluded a posteriori by identifying those cases where ¢ -= u c 0 changes sign (no critical level is considered here).
The eigenvalues associated to modes A and B are real and correspond to realistic values of phase speed for propagating waves, although in many cases the eigenvalue problem tends to overestimate negative phase speeds (this is thought to be a consequence of barotropicity). However, in a small number of cases, among the solutions of (12), couples of complex eigenvalues were found. It was observed that the number of such solutions decreases as the input zonal wind profile becomes smoother and increases as its strength grows. This class of solutions, interpretable as related to barotropic instability, were not considered being outside the scope of this study.
The meridional structure f ( ) p of the second order solution is obtained by solving the differential equation (14) with a finite difference scheme. The nonlinear self-interaction coefficient χ is finally computed by numerical integration of (17).
Individual wavenumber components are obtained through Fourier analysis by calculating the complex functions  was kept constant and f 0 was varied in the range 15°N-35°N for k=1-3 and for all meridional modes. Not unexpectedly, variations of the latitudinal domain extension proved to affect the estimate of the eigenvalues c. Recalling the classical expression for the phase speed of free-propagating barotropic Rossby waves, it appears plausible that the values of c decrease along with the extension of the latitudinal domain. As the domain narrows, the meridional wave number increases and the phase speed tends to approach zero. It is also worth noting that, as the meridional structure of the asymmetric field becomes more complex, its phase speed decreases.
According to the theoretical predictions of the adopted multiple-scale analysis approach, the general solution of the reduced barotropic equation at the second order in the perturbative parameter ε is given by (20) and is characterized by the presence of a correction to phase speed given by ε 2 χ. The signature of self-nonlinearity on the propagation of planetary scale waves was sought in ERA-Interim data by looking for a positive correlation between measured wave phase speeds and the values of ε 2 χ obtained from theory, with ε given by (25). However, it can be hardly expected to find statistical evidence of such a correlation in reanalysis of data, since eigenvalues and eigenfunctions depend on the details of the zonal wind profile, which change from day to day. In the scatter plots of the measured phase speed c * versus computed nonlinear corrections ε 2 χ, a peculiar feature is observed: as manifest in figure 2, large (positive and negative) values of the nonlinear corrections to the phase speed tend to occur in concomitance with very small or zero c * for modes A and B. Also evident is the accumulation of fast traveling waves on the axis of zero self-interaction. This result supports the hypothesis that self-nonlinearity preferentially acts with the effect of slowing down wave propagation, and that its role is fundamental especially in the dynamics of stationary (or nearly stationary) waves.

Conclusions
A multiple-scale analysis of the solution for the barotropic vorticity equation was performed in absence of any external forcing or topography, deriving a nonlinear self-interaction second order correction to phase speed proportional to the square of wave amplitude.
The above formulation was used for analysing zonal wind data and geopotential fields at the 500 hPa pressure level from ECMWF ERA-Interim reanalysis in the latitudinal domain 25.5°N-90°N. The analysis focused on months ranging from October to April in the period 1979-2015, as during these months ultra-long fluctuations activity is prominent in the Northern Hemisphere. The estimated statistics of nonlinear self-interaction shows that its strength tends to increase with the complexity of the structure of the fluctuations, depending on their wavenumber and meridional profile. Although the estimated corrections to phase speed are often small or negligible, results show that higher values of self-interaction are found for stationary or quasi-stationary fluctuations, supporting the hypothesis that this nonlinear process plays an important role in the dynamics of the standing part of the Hayashi spectra.
The use of the barotropic vorticity equation, although traditionally justified by the equivalent barotropic nature of large-scale waves, is a limitation especially for what concerns the correct estimation of the theoretical phase speed of Rossby waves. Moreover, it is well known that the energetics of planetary waves is dominated by baroclinic conversion [30], i.e., the conversion from zonal to eddy available potential energy associated with a small vertical tilt of the planetary waves. The incorporation of these features in the theory is possible and can be tackled by extending the perturbative approach to the framework of a baroclinic atmosphere.
A final, important aspect implicit in this study is the problem of the meridional confinement of planetary waves, which is a physical requirement at the basis of the theory presented here. Time fluctuations of the zonal wind (here neglected) have the statistical effect of confining the propagation of Rossby waves [31], and constitute an important feature to be taken into account in future generalizations of the self-interaction theory. a Figure 2. Observed phase speed c * versus nonlinear phase speed correction computed for modes A (a) and B (b). Results for wavenumbers k=1, 2, 3 are shown together. Units are dimensionless (one unit corresponds to 10 ms −1 ).