Dynamics

Abstract. It is paradoxical that, while atmospheric dynamics are highly nonlinear and turbulent, atmospheric waves are commonly modelled by linear or weakly nonlinear theories. We postulate that the laws governing atmospheric waves are in fact high-Reynolds-number (Re), emergent laws so that – in common with the emergent high-Re turbulent laws – they are also constrained by scaling symmetries. We propose an effective turbulence–wave propagator which corresponds to a fractional and anisotropic extension of the classical wave equation propagator, with dispersion relations similar to those of inertial gravity waves (and Kelvin waves) yet with an anomalous (fractional) order Hwav/2. Using geostationary IR radiances, we estimate the parameters, finding that Hwav a 0.17 ± 0.04 (the classical value = 2).


Introduction
The atmosphere is a highly turbulent system with the ratio of nonlinear to linear viscous terms -the Reynold's number -typically of the order ≈10 12 . At the same time, there is no doubt that atmospheric waves exist and play an important role in transferring energy and momentum. These empirical facts only become problematic when we consider the numerous apparently successful studies comparing data with linear (or weakly nonlinear) theory, commonly (for gravity waves) with the Taylor-Goldstein equations or with the linearized shallow water equations. For example, in the words of (Nappo, 2002) "Almost all of what we know about the nature of gravity waves is derived from the linear theory" (emphasis in the original).
Although one may easily get the impression that linear wave theories have been empirically confirmed, a closer look reveals that what has typically been scrutinized are the linear theory dispersion relations. Considering the example of gravity waves, we find that these have mostly been tested in the horizontal (and occasionally in the vertical) directions. Other predictions of the corresponding linear theory -"polarization relations" -are invoked but are only used in a diagnostic mode so that they cannot be considered to have been convincingly validated. Recently, linear gravity wave theory has been directly brought into question by data from dropsonde pairs. For example, such pairs have directly shown that certain terms neglected in the Taylor-Goldstein equations are typically much larger than the corresponding retained terms . Also the pairs have clearly shown that the vertical structure of the atmosphere is composed of a fractal hierarchy of unstable layers through which linear gravity waves would not able to propagate (Lovejoy et al., 2008a).
One argument for justifying the application of linear wave theories to the atmosphere are the apparent existence of large laminar flow regions. However, since the 1980's -and largely thanks to the development of multifractal cascade modelsthere has been dramatic progress in understanding atmospheric intermittency. It is now clear that a prime characteristic of fully developed turbulence is that most of the important fluxes are concentrated in highly sparse (fractal) sets so that much of the flow appears relatively calm. The modern understanding is that by its very nature, turbulence is highly intermittent so that on any realization of a turbulent process, there will be violent regions in proximity to ones of relative calm. However, examination of the apparently calm regions shows that they also have embedded regions of high activity and as we zoom into smaller and smaller regions this strong heterogeneity continues in a scaling manner until we reach the dissipation scale. This explains why aircraft measurements of the wind invariably find roughly k -5/3 (i.e. turbulent) spectra even in apparently calm regions. Large scale regions of true laminar flow have yet to be documented by actual measurements. On the contrary, the multifractal, multiplicative cascade picture has been well verified even at large scales (e.g. Lovejoy and Schertzer, 2010). Therefore, it would be a mistake to separate these regions of high and low "turbulent intensities" and associate them with different mechanisms or to apply nonturbulent (linear) wave theories to regions of apparent calm.
In the last few years, (nonlinear) scaling theories of waves have become more compelling. This is because empirical evidence and theoretical arguments have amassed to the effect that atmospheric dynamics give rise to emergent high Reynolds number scaling laws with different horizontal and vertical exponents. This allows the horizontal scaling to accurately apply over huge ranges in scale, (see Lovejoy and Schertzer, (2010), Lovejoy and Schertzer, (2013) for recent reviews). Although based on the classical laws of turbulence, they involve extensions to account for (multifractal) intermittency and anisotropy. Their success underlines the fundamental role of scale symmetries in constraining the high Re dynamics. All this motivates the question: are atmospheric waves also scaling turbulent phenomena? If this is the case, we may logically expect anomalous wave propagators that could readily have dispersion relations identical or nearly indistinguishable from their classical counterparts, while simultaneously having nontrivial consequences for the dynamics and for our understanding.
If dispersion relations from linear theory and those from strongly nonlinear theory can be very similar to each other, then how might one empirically distinguish them? The obvious way is to note that linear theory also predicts the entire space-time propagators relating the wave forcings and responses. A key characteristic of linear theories is that they involve integer powers of the (space and time) differential operators, and this strongly constrains the form of the propagators; below, we show how this allows us to test the theory by seeking possible anomalous propagator exponents. Below, we investigate this using geostationary satellite infra red radiances.
The following paper attempts to show how scaling propagators with both turbulent and wavelike characteristics could arise while being consistent with both (anisotropic) turbulence theory and observations. However, let the reader be warned that while the turbulent part of the propagator -which was derived and empirically tested elsewhere (it is summarized here in an appendix) -is reasonably well grounded, in contrast, the wavelike part -the subject of this paper -is fairly speculative, it is perhaps little more than a proof of concept. On the theoretical side, the reason is that with only scaling symmetries to guide us, the possibilities are very broad while on the empirical side, over the scaling range (120-5000 km in space and 3-100 hrs in time) the turbulent part of the spectrum is by far the dominant one accounting for an empirical range of spectral densities of ≈ 10 5 leaving the (residual) wavelike part to account for the remaining factor of 2-4 (the mean residual is 0.9±0.5) in the dynamical spectral scaling range.

Fractional propagators and turbulence
In order to motivate our model, let's consider: the classical wave equation for the wave I with forcing f: V is the wave velocity, r is the position vector and t the time variable.
As usual, we can solve eq. 1 by taking Fourier transforms (denoted by tildas): where k is the wavevector, ω the frequency and ( , ) g k ω is the propagator. This propagator is symmetric with respect to an isotropic space-time scale transformation by factor λ −1 : However we anticipate that at high Re "effective propagators" may emerge constrained by the same overall scaling symmetry but with some other "anomalous" exponent H≠2. In this case we obtain fractional propagators corresponding to fractional generalizations of the wave equation: Although we will only require fractional propagators, if needed, we could define the fractional differential operator in eq. 4 by the inverse Fourier Transform of g (k,ω ) −1 (or see e.g. (Miller and Ross, 1993) for fractional differential equations). If we seek the real space solution of eq. 1 or 4, we can use the fact that Fourier space products (eq. 2) correspond to real space convolutions ("*") hence the solutions to eqs. 1, 4 are: I (r,t) = g(r,t) * f (r,t) so that the propagator links the forcing (f) to the response (I).
In order to estimate the g(r,t) we can appeal to the method of stationary phase (e.g. Bleistein and Handelsman, 1986) which ensures us that the dominant contribution to g(r,t) is due to the wavenumber-frequency region over which g (k,ω ) is singular, this singularity defines the dispersion relation and accounts for its origin and significance. For both the classical eq. 1 and the nonclassical eq. 4, we find the dispersion relation: which is therefore of fundamental importance, a fact which is true for any H>0, not only for positive even integer values of H.
Before attempting to estimate propagators of real data, we must take into account the fact that atmospheric waves occur in the presence of turbulence. Indeed the spectrum is so strongly dominated by a "turbulent background" that it must first be removed before evidence of any wavelike propagator can be observed. This is paradoxical since the wavelike part implies the existence of a singular set over a surface in (k x , k y , ω) space and should be easy to detect (actually, the topology need not be so simple, see section 5). However, the singularity is of sufficiently low order and spectral estimates are sufficiently noisy that in practice the surface is hard to observe. Indeed it is much easier to study 2-D subspaces obtained by integrating out one of the spectral coordinates (which also reduces the "noise"), although if H is small enough this can integrate out the singularities. Indeed, one of the main techniques for empirically investigating atmospheric waves ( (Wheeler and Kiladis, 1999 ;Hendon and Wheeler, 2008 ;Dias et al., 2012)) integrated over k y space to yield a (k x , ω) 2-D space spectrum while also using an ad hoc averaging technique to removing the turbulent.
Following (Wheeler and Kiladis, 1999), we also use IR imagery (but at hourly not daily resolution), we use instead a theoretically motivated turbulent spectrum to search for evidence of anomalous wave propagators. To understand this, recall the classical Kolmogorov law of three dimensional isotropic turbulence: where I is a component of the wind, ΔI is a fluctuation, Δr is a vector displacement, ε is the turbulent energy flux and the equality is understood in a statistical sense. In Fourier space this becomes: comparing this with eq. 2, we see that ϕ k ( ) is the forcing and g tur k ( ) is the spatial part of a (fractional order) propagator (a Green's function). Now recall that for real I, , if in addition we assume statistical translational invariance ("statistical homogeneity"), then we may define the spectral densities P I , P ϕ , by: where "< . >" denotes ensemble averaging.
To obtain the classical Kolmogorov-Obhukhov k -5/3 law we use: where P 0 is a dimensional constant, d is the dimension of space and K(2) is the second order intermittency correction. This yields: The angle integrated ("isotropic") spectral density E(k) is then given by integrating P over shells in Fourier space. Ignoring constant factors (2π in d = 2, 4 π in d = 3), we obtain the (intermittency corrected) isotropic Kolmogorov law: (since H=1/3, we see that the nonintermittent K(2)=0 law does indeed have exponent β=5/3).
A basic consequence of wide range spatial scaling of atmospheric fields (in particular the wind) is that the spectrum and spectral density of the turbulent fluctuations in horizontal wavenumber -frequency (k x , k y , ω) space follow the straightforward space-time extension of eq. 6: where P I (k,ω ), P ϕ (k,ω ) are space-time spectral densities, ( , ) tur g k ω is the turbulent propagator. To obtain the form of ) , ( ω k g tur , we follow (Lovejoy and Schertzer, 2010 ;Pinel et al., 2013) outlined in appendix A (see eq. A13) to obtain: where: where k is the (horizontal) spatial (Fourier) scale function and a is the north-south/ east-west aspect ratio; when a=1, we obtain k k = (note that if needed, more complex spatial scale functions may be used, these replace the vector norms in the isotropic theories). The transformation ω ω′ → combines the effects of a mean advection by velocity μ and the statistical variability of the advection wind about its mean is accounted for by σ. Note that a) the factor i in eq. 12 is necessary so that the propagator respects causality, and b) overall tur g respects the same isotropic scaling symmetry as the wave propagator, eq. 3 but with exponent H tur .

Fractional propagators and waves:
With the exception of the weak singularities associated with waves, the turbulence dominates the spectral density, the P I (k,ω ) given in eq. 11 with the propagator eq. 12 already gives a good approximation to the empirical spectral density. This may be seen in fig. 1 using MTSAT data (described below) which shows the 1-D spectral densities E(k x ), E(k y ), E(ω) obtained by successively integrating out various pairs of variables from ( , , ) Pinel et al., 2013). This log-log linearity on this figure directly shows that the spectra are scaling and near perfect superposition of the 1-D spectra shows that the scaling exponents are essentially identical so that (in conformity with the conclusions of appendix A and the form eq. 12) the radiance field structure functions are symmetric with respect to isotropic scale This turbulence part corresponds to the "background" spectrum obtained by Wheeler and Kiladis, (1999); any wave behaviour is to be found in deviations from this.
A simple model that takes into account waves while respecting both the spacetime scaling and the turbulent forcing and background is obtained by including a factor in the overall propagator. To be "wave-like", wav g wav g must be causal, unlocalized in space-time and must also be chosen so that the overall scaling symmetry of the system (eq. 3) is respected by ( ) , I g k ω . Following Wheeler and Kiladis, (1999) who factored the spectral density into a "red noise" turbulent background and a wave part and inspired by eq. 2 and eq. 12, we can use the form: with given by eq. 12 and given by: this is a generalization of eq. 2 to account for spatial anisotropy (with k k → ). The replacement k ω ω → + ⋅μ is the classical advection transformation (see e.g. Nappo, 2002); as in the turbulent propagator, we have included the extra factor σ to take into account the statistical variation of the advection velocity. Finally, the parameter v wav is the phase speed nondimensionalized by the turbulent velocity V w (eq. 12). Note that the overall propagator satisfies the scaling symmetry eq.3 with H = H I g tur +H wav .
Due to , the overall propagator wav g I g yields the dispersion relation: With respect to the "background" advection (μ), wav v σ is the effective wave speed which takes into account the true mean wave speed (v max ) and the statistical variability via σ. By taking appropriate scale functions k one can obtain dispersion relations close to gravity and other waves (see Lovejoy et al., 2008b).
Of more relevance here are Kelvin waves which are the low Coriolis parameter/high "effective thickness" limit of the inertial gravity (Poincaré) wave dispersion relations often invoked at these space-time scales. First, for only one spatial (zonal) dimension, we may note that Kelvin waves are a special case of eq. 17 Considering the full horizontal plane, Kelvin waves are "channelled", only propagating in the zonal direction. To obtain some channeling while maintaining the same overall scaling symmetry, we could replace use the spatial (Fourier) scale function ( ) In eq. 18, we have followed the assumption in the isotropic case (eq. 9) that the forcing of the flux has the same scale symmetries as 2 tur g ; from eq. A13 we see that is the spectral exponent of the flux and, P 0 a dimensional constant determined by the climatological (low frequency) average forcing.

Data Analysis
We follow (Wheeler andKiladis, 1999, Hendon andWheeler, 2008)  We now consider the three 2D spectra, obtained by successively integrating the 3D spectral density over k x , k y and ω. The fit is sufficiently good that we can use the above regression with H wav = 0 to estimate all the turbulent parameters. However for the 1-D spectra to have fixed exponents, when fitting the wave part we must use the constraint H = H tur +H wav so that the 1-D spectral slopes are not affected. In this way we find an optimum relative weighting for the turbulence and wave contributions. in the three 2D spaces (excluding the diurnal spikes and the origin) which is small considering that the 2D space signal P(k, ω) varies over ~ four orders of magnitude.
The orientations of the contours of P(k x , ω), P(k y , ω) is a consequence of the non zero In order to directly visualize the wave behaviour, Wheeler and Kiladis, (1999) removed a "turbulent background" (estimated with an ad hoc averaging technique) from their (k x , ω) 2-D spectrum and tried to identify maxima in the residual with linear theory dispersion relations. Inspired by their work, we removed (by dividing) from the empirical 3D spectral density the "turbulent background" estimated from the fit of eq.
18 with H wav =0 (i.e. the purely turbulent spectral density from which we obtained figure 2a). The residual from which wave behaviour is to be identified (and which is to be described by the wave part of eq.18) is presented in figure 3 for the (k x , ω) 2-D space with ω>0 (i.e. after integrating over k y ). We observe a region of maxima ( fig.3 in grey) for k x >0 which is similar to the residual obtained by Wheeler and Kiladis, A key point about Kelvin waves and fig.3 is that they are asymmetrical in the zonal direction (k x ). In contrast, the simple form of wav g used in eq.18 involves a Fourier space scale function ( 1/ 2 2 2 2 x y k k a k = + ) which is symmetrical in k and which involves, in the (k x , ω) space, maxima lines (coming from the singularities) for k x >0 as well as for k x <0, which is incompatible with figure 3 (grey region). However, the only constraints on the form of are that it must respect causality, that g is real, (hence Replacing in eq.18 by eq.19 preserves the quality of the fit of the total (turbulent-wave) spectral density, (see the 2D subspaces in figure 2c) and gives a spectra wav g ( , ) wave x P k ω close to the data, including a maxima line which is close to the maxima in the residual presented in figure 3.

Refined singularity analysis
The above analysis is paradoxical since our hypothesis is that there is a singular surface in (k x , k y , ω) space yet analysis of the 1-D and 2-D sections showed no direct evidence of singular behaviour. This is consistent with the finding that 0<H wav <1 implying that the singularities are integrated out in the lower dimensional sections. In order to display potential singularities, we are therefore forced to study the full 3D density P(k x , k y , ω) recognizing that most of the variation is due to the turbulent part and that the wave part -being only weakly singular -is expected to manifest itself in local maxima, perhaps with surface-like topology (line-like in 2D sections). To this end, we implemented an ad hoc singularity detection algorithm that "scans" parallel to the axes to estimate local maxima successively in the k x , k y directions. In principle considering the maxima in a single direction is adequate, but in practice the singular surface has parts that are roughly parallel to a given axis; the resulting ambiguity can be resolved by determining maxima in two orthogonal directions.
The results are shown in fig. 4, where we compare such an analysis with the theoretical behaviour (with from eq.19) for constant ω sections. A drawback of the method is that it does not distinguish maxima due to the turbulent contribution and from the (presumed) wave contribution and in the empirical case, the separation is not always evident. In the figure, the two have been distinguished by the color of the lines. We see that although far from perfect, the semi ellipses indicating the wav g theoretical singularity (dispersion relation) are close to the empirical ones. Given that we used a straightforward generalization of the classical wave equation with only one new parameter v wav (two if we include H wav = 0.08, but this doesn't affect the singular surface). Given that the wave part of the overall propagator eq. 16 as well as the postulated multiplicative decomposition (eq. 15) are not much more than the simplest analytical hypotheses, the results are quite encouraging, yet they indicate some of the difficulties.

Conclusion
The atmosphere is highly nonlinear yet displays both turbulent and wavelike behaviour over huge ranges of space-time scales. Theories explaining the turbulent aspects assume that the dynamics are strongly nonlinear and scaling, in contrast, the corresponding wave theories have been linear or weakly nonlinear. We proposed that the paradox can be explained by noting that although linear theory predicts propagators, only the relations implied by the singular part of the latter -the dispersion relations -have been tested to any extent. However, linear theories invariably involve integer differential operators and corresponding integer ordered propagators so we may test these theories by examining the propagators -or at least their squared moduli which are amenable to empirical spectral determination.
The mathematical structure of the turbulent laws that link the observables to driving turbulent fluxes (such as energy fluxes) use scaling (turbulent) propagators which are very similar to that of wave equations except that the latter are singular. To account for both waves and turbulence, the actual propagators need only respect scale symmetries and can be modelled as products of turbulent-like and wave-like (spacetime localized and unlocalized) propagators with both involving anomalous exponents.
The wave propagator we used involves the mean horizontal turbulent wind and energy flux as well as a mean background wave advection velocity, obtained as an (anisotropic, fractional) generalization of the classical wave equation (which is approximately satisfied by inertial gravity waves and Kelvin waves). Using two months of MTSAT hourly geostationary IR radiances, using a simple form for the wave propagator, we found that the best fit involved an anomalous wave scaling exponent H wav ≈ 0.17±0.04; for comparison, the classical wave equation has the integer value H wav = 2.
Investigating the wave structure, we followed Wheeler and Kiladis, (1999), dividing our empirical spectral density by a theoretically estimated "turbulent background"; the maxima in the residual are to be identified with wave dispersion relations. Since there are very few theoretical constraints on the form of the wave propagator, we chose a simple ansatz that is compatible with the observations. This paper is simply an early attempt to understand waves in highly turbulent media using scaling symmetries as constraints. On the one hand, these symmetries are so broad that they provide only limited guidance, on the other hand, the turbulent part -with no wave contribution at all -explains almost all of the observed dynamical spectral scaling range of factor ≈ 10 5 leaving only a small residual (a factor 2-4) for the wavelike part. The empirical situation would certainly be improved if data from other fields covering wide scale ranges in the full (x, y, z, t) space could be found, but at present the (x, y, t) geostationary radiances are apparently the best available.
Therefore, this paper should be seen more as a proof of concept than as providing definitive results. The main conclusion is thus that strongly turbulent atmospheric dynamics are a priori compatible with the observed waves, that to understand them, that one doesn't need to invoke the existence of large laminar regimes and linear theories.

Appendix A: The space-time turbulent spectrum
The 23/9D model of spatial turbulence Schertzer and Lovejoy, (1985a,b) involves wide range scaling in the horizontal and vertical directions but with different scaling exponents; the horizontal being dominated by energy fluxes, and the vertical by buoyancy variance fluxes. Since the infra red radiances are essentially (x, y, t) (horizontal-time ) fields we needn't explicitly consider the vertical, however we do need to extend the model to space and time. In this appendix, we summarize the arguments developed in Lovejoy et al., (2008b), Schertzer, (2010, 2013) and in Pinel et al., (2013).
The first step is to rewrite the isotropic eq. 5 in a more general anistropic scaling manner by replacing the vector norm by a space-time scale function R Δ : ( ) where we have used the subscript R Δ to emphasize that the flux is at resolution R Δ . The space-time scale function is symmetric with respect to generalized scale changes T λ : where G is the generator of the scale changing operator T λ . Ignoring for the moment advection, and taking I as a horizontal velocity component then the canonical (simplest) nondimensional scale function compatible with the Kolmogorov law is: where H t = (1/3)/(1/2) = 2/3 and L w , τ w are the outer scales of the scaling in space and in time and a is a north-south/east-west aspect ratio. The outer scales are linked by the overall average energy flux ε: τ w = ε -1/3 L 2/3 (τ w is the lifetime/"eddy-turn-over time" of structures size L w ). Successively substituting into eq. A3 and the latter into A1, yields the Kolmogorov law in the horizontal directions (e.g. ) and time, the Lagrangian law (Δv = The next step is to consider the scale function corresponding to a constant advection in the horizontal Due to Gallilean invariance, (i.e. under ; ; x y Since there is no scale separation, v is a turbulent velocity, so that over a given region, it will be dominated by the advection due to the largest eddies. For the same reasons, the link between L w and τ w is via a turbulent velocity, a consequence of which is that the pure temporal development (Lagrangian) term may be neglected.
The "effective" scale function we seek is thus obtained by averaging over the turbulence (this argument is not completely rigorous since due to the intermittency, averaging over other powers of R Δ will give somewhat different parameters). Turbulence will have two effects on eq. A5: the mean advection by v, and the effect of turbulent variability. The overall result (for more details, see Pinel and Lovejoy, 2013) is: where the "effective" scale function is given by: where we have used the nondimensional variables: is the overall mean advection over the region considered and is the large-scale turbulent velocity at planetary scale.
The intermittency corrections come from the scaling of the flux ϕ: so that overall: ( ) The structure function exponent ξ(2) thus takes into account the scaling exponent H as well as the intermittency correction.
To obtain the corresponding spectral density ( , ) I P k ω (and hence, we follow the development presented by Lovejoy and Schertzer, (2010) (see also Pinel et al., (2013)) and use the general relation between structure functions and spectra: so that the effective real space scale function defines an effective Fourier space scale function: With this scale function, we have the nondimensional spectra: (with d=3 for horizontal space-time).
(2) (2) -1 with wav g from eq.19); the blue, the empirically estimated curve using the ad hoc algorithm and the green and red show maxima but presumed to originate in the turbulence "background" (they are very close to the axes).