On the nonlinear shaping mechanism for gravity wave spectrum in the atmosphere

The nonlinear mechanism of shaping of a high vertical wave number spectral tail in the field of a few discrete internal gravity waves in the atmosphere is studied in this paper. The effects of advection of fluid parcels by interacting gravity waves are taken strictly into account by calculating wave field in Lagrangian variables, and performing a variable transformation from Lagrangian to Eulerian frame. The vertical profiles and vertical wave number spectra of the Eulerian displacement field are obtained for both the case of resonant and non-resonant wave-wave interactions. The evolution of these spectra with growing parameter of nonlinearity of the internal wave field is studied and compared to that of a broad band spectrum of gravity waves with randomly independent amplitudes and phases. The calculated vertical wave number spectra of the vertical displacements or relative temperature fluctuations are found to be consistent with the observed spectra in the middle atmosphere.


Introduction
Observations of the mesoscale wind velocity and temperature fluctuations in the middle and upper atmosphere often reveal a presence of a small number of discrete gravity waves (Sica, 1999;Gurvich and Chunchuzov, 2005).Other studies suggest that the fluctuations are induced by a broad band spectrum of gravity waves (Allen and Vincent, 1995;Hines, 1991).To explain the observed forms of the vertical wave number spectra of temperature fluctuations in the stably stratified layers of the atmosphere (troposphere, stratosphere and Correspondence to: I. P. Chunchuzov (igor.chunchuzov@gmail.com)mesosphere) we analyze in this paper a nonlinear shaping mechanism for these fluctuations.
For the ocean a similar shaping mechanism was proposed earlier by Allen and Joseph, 1989 (hereinafter AJ89), who tried to find a physical explanation for the well-known empiric forms of the oceanic internal wave spectra derived by Garrett and Munk (1975).This mechanism was associated with the influence of the wave-induced advection of fluid parcels on both spatial and temporal internal wave spectra.The AJ89' approach was later applied by Chunchuzov (1996) to the atmosphere, and then developed by Chunchuzov (2001), Hines (2001) (hereinafter HO1), and Chunchuzov (2002) (hereinafter CO2) with substantial modifications of the AJ89'approach (the cause of these modifications and their detailed description are given in HO1 and CO2).Based on the nonlinear shaping mechanism we study here the dependence of the spectral forms on whether wave sources excite a few discrete gravity waves during observational period, or a broad band spectrum of random waves.
A nature of the mesoscale fluctuations in the atmosphere is debated for a long time (see, for instance, Fritts and Alexander, 2003, Sect. 4.1).In certain range of high vertical wave numbers k z the observed vertical wave number spectra of the temperature fluctuations show a k −3 z -power law.It was noticed by Hines (1991) that the gravity waves with the high k z , which are within the k −3 z -spectral tail, have very low horizontal phase speeds, comparable to the horizontal wind velocity fluctuations induced by the waves themselves.Such waves should strongly interact due to advective nonlinearity of the Eulerian fluid motion equations.
To take into account the wave-induced advection of fluid parcels it was suggested in a number of works to use a Lagrangian frame of variables for studying the dynamics of gravity waves, and then perform a variable transformation to the Eulerian frame, where the advection takes place (AJ89; Chunchuzov, 1996Chunchuzov, , 2001;;Eckermann, 1999;HO1;CO2;Broutman et al., 2004;Pinkel, 2008).As an exact Published by Copernicus Publications on behalf of the European Geosciences Union.
transformation it strictly takes into account the advective effects associated with the nonlinear terms (v∇)v in the Eulerian equations of motion without using any approximations for these terms.Using a Lagrangian approach it was found in CO2 that a strong nonlinearity of the wave field generates a 3-D Eulerian wave number spectrum with a k −5 -power law decay at high wave numbers k.This spectrum is of highly anisotropic form as a result of a balance between the nonlinear wave energy transfer from the characteristic (vertical and horizontal) scales of internal wave sources toward smaller vertical and larger horizontal scales, and the dissipation of wave energy at small vertical scales due to wave breaking processes.Such cascade-like energy transfer in the 3-D wave number space is caused by non-resonant wave-wave interactions, which along with wave energy dissipation play a key role in shaping of the equilibrium gravity wave spectrum.The wave-like fluid motions generated by the non-resonant interactions resemble anisotropic and vertically oriented vortices rather than the linear gravity waves, because the dispersion surfaces of these waves in the frequency-wave number space are completely "smeared" by advection.
The hypothesis about a significant role played by forward energy cascade in shaping of the spectra of the mesoscale fluctuations in the atmosphere was earlier proposed by Dewan (1997) on the basis of his saturated-cascade similitude theory.Using this theory Dewan found the forms of the horizontal and temporal gravity wave spectra, which were close to their observed forms, although he traditionally interpreted the vertical spectra as a result of the saturation of linear gravity waves caused by their convective or shear instabilities.The merging of wave saturation and wave cascade processes was assumed to be possible due to the existence of a unique relation between the vertical wave length and period of gravity waves.Lindborg (2006) also used a similitude theory added by numeric simulations of energy cascade in stratified fluid.This allowed him to derive the k −3 z and k −5/3 h forms for the vertical spectra (scales from 100 to 1000 m) and the horizontal spectra (scales from about 1 to 500 km), respectively.However, contrary to the linear wave saturation hypothesis, Lindborg (2006) assumed that both the horizontal and vertical spectra arise ". . .from one and the same type of nonlinear chaotic motion. . ." governed by the fully nonlinear Boussinesq equations.Such highly anisotropic motions significantly differ from those induced by linear gravity waves, therefore he suggested for these motions a neutral term "layer" instead of "wave".The conclusion made by Lindborg (2006) about the importance of a nonlinear energy cascade in shaping of highly anisotropic mesoscale fluctuations in stably stratified fluid is in agreement with that previously obtained in CO2 based on Lagrangian approach.
Recently, the 3-D spectrum obtained in CO2 was slightly modified by Gurvich and Chunchuzov (2008) who assumed that the anisotropy of the temperature fluctuations decreases with a decrease of their vertical scale.Based on this assump-tion they developed the model of the 3-D spectrum of temperature fluctuations, from which both the vertical and horizontal wave number spectra were obtained.These spectra were shown to be consistent with their observed forms in the middle atmosphere.
The nonlinear shaping mechanism proposed in CO2 forms a certain Eulerian frequency spectrum of the temperature and wind velocity fluctuations at a fixed point of space (Chunchuzov et al., 2006).This spectrum was derived from the 4-D frequency-wave number spectrum obtained in CO2.It was found that at high frequencies (ω>N ), which are beyond the buoyancy frequency N, the nonlinear wave-wave interactions generate in the frequency spectrum the ω −3 -tail of the same nature as the k −3 z -tail in the vertical spectrum.For the intermediate range of frequencies (f ω<N,f is inertial frequency) the frequency wind velocity spectrum was shown to follow the εω −2 power law, where ε=σ 2 N/χ is the mean generation rate of wave energy from random gravity wave sources, σ is the rms value of the wave-induced horizontal wind velocity fluctuations, χ=σ/σ w characterizes the anisotropy of the Lagrangian wave field, σ w is the rms value of the vertical wind velocity fluctuations.In this range the frequency spectra of the wind velocity fluctuations measured in stably stratified lower troposphere (by acoustic anemometers and acoustic sounding) in average (over all the measured spectra) decayed with growing ω as εω −2 (Chunchuzov et al., 2006).
Recently, the Lagrangian approach, but slightly simplified, was applied to the oceanic waves by Pinkel, (2008) to explain the observed velocity shear spectrum in the ocean.He considered the time-varying advection of fluid parcels caused by a few energy contained sinusoidal constituents of the wave field such as inertial and tidal components.The latter were taken in so-called semi-Lagrangian frame.The kinematic distortion of the wave field due to transfer from this frame to the frame of measurement platform resulted in the advective "smearing" of the discrete spectral lines.Analyzing this effect Pinkel (2008) assumed that the apparently continuous wave number-frequency spectrum of the oceanic shear can result from a wave field populated at only a few intrinsic frequencies.It was noticed earlier in CO2 that the increasing with wave amplitudes kinematic effects of advection become significant only when the nonlinearity of the wave field in the Lagrangian frame becomes significant as well.Therefore, the shaping mechanism for the wave spectrum considered here is not associated with a purely kinematic effect, but results from a combined effect of both the advective and dynamic nonlinearities of the wave field.
The important result of the nonlinear theory developed in CO2 was the obtained form for the 4-D gravity wave spectrum, which allowed one to explain the observed forms of both the frequency and spatial 1-D (vertical and horizontal) spectra based on the same shaping mechanism.The generality of this mechanism was shown by parameterization of the vertical, horizontal, frequency spectra and the corresponding second-order structure functions (under nonlinear saturation of the fluctuations) through the same main parameters N,σ and χ, and by the agreement of the predicted and observed spectra (Gurvich and Chunchuzov, 2008;Chunchuzov et al., 2006).However, the forms of other statistical characteristics of the fluctuations such as their probability density functions and third-order structure functions still need to be obtained from the same theory and compared to the experimental data available.
The nonlinear theories of wave spectra proposed by AJ89, Chucnchuzov (1996), HO1, Hines (2002), and CO2 were recently tested by Klaassen and Sonmor (2006) (hereinafter, KS06) and Klaassen (2009), who came to the conclusion that their kinematic model does not provide any support for these theories.KS06 found that for the typical wave amplitudes required to produce the k −3 z -spectral tail the vertical profiles of the Eulerian displacement wave field become multivalued and the transformation from Lagrangian to Eulerian variables becomes invalid.As an alternative to the nonlinear shaping mechanism for this tail, Klaassen (2009) suggested that ". . .through instability, atmospheric wave fields create and possibly coexist with a field of smaller-scale (perhaps turbulent) eddies, which act to exert a force on the mean background flow.In other words, the secondary field of eddies acts as a dissipation mechanism for the larger scale internal waves, producing momentum deposition and perhaps saturation as the latter continue to propagate. . ." It is necessary to note that the multivalues can arise in the profile of the wave field, because the latter is a solution of nonviscous fluid motion equations.Another well-known example of such solutions is a plane acoustic wave of finite amplitude, whose wave profile distorts with increasing distance from a plane source until it becomes multivalued near its wave front (Lighthill, 1978;Rudenko and Soluyan, 1977).By taking into account a molecular viscosity and thermal conductivity one can stabilize the nonlinear steepening of the wave profile and prevent the arising of multivalues.In the case of internal waves propagating through realistic atmosphere the different types of wave-induced instabilities may prevent the arising of discontinuities in the wave field profile through wave breaking processes and transferring wave energy into the energy of turbulent eddies.However, these processes as shown in CO2 occur only within the thin spatial regions whose thickness is small compared to the vertical scales, 2π/k z , typical for the k −3 z -tail.In these local regions of space a sink of wave energy balances the nonlinear wave energy transfer through the entire spectral tail.
Since wave drag parameterization schemes are sensitive to the height dependence of the forms of the vertical wave number spectra and of the characteristic scales of wave breaking processes the estimating of these scales is of great importance for solving the problem of parameterization of gravity wave drag in the atmospheric circulation models (McLandress, 1997;Fritts and Alexander, 2003).Beside this, the knowledge of the space-time spectrum of the wave-induced temperature and wind speed fluctuations is needed for the prediction of the statistics of the amplitude and phase fluctuations of low-frequency acoustic waves (Chunchuzov et al., 2006;Ostashev et al., 2005) and electromagnetic waves (Gurvich andChunchuzov, 2003, 2005), propagating through a realistic atmosphere.
In this paper we consider the same nonlinear shaping mechanism for gravity wave spectrum as in CO2, but for only two waves.Based on the Lagrangian fluid motion equations, derived in Sect.2, we will study in Sect. 3 a nonlinear generation of the harmonics of two given waves in the Lagrangian frame of variables, and calculate their displacement field by using perturbation method.In Sect.4, the obtained field and its vertical wave number spectrum will be studied in the Eulerian frame.Both resonant and non-resonant wave-wave interactions will be analyzed.Such study will be extended in Sect. 5 to the case of only one ducted gravity wave excited by wave sources to show a generality of the shaping mechanism discussed here.In Sect.6 the spectra for the case of a few discrete waves will be compared to those for the case of a high number (n 1) of gravity waves with randomly independent amplitudes and phases.The new forms obtained for 3-D and 1-D wave number spectra will be analyzed.

Equations for nonlinear internal wave field in the Lagrangian frame of variables
We start with a system of motion equations for nonviscous, stably stratified, and nonrotating fluid in Lagrangian frame of variables r=(a,b,c) and t (Lamb, 1932;Gossard and Hooke, 1975): where g is gravity acceleration, γ is adiabatic constant, x,y,z are the co-ordinates at a moment t of a fluid parcel, whose undisturbed position under static state of the atmosphere is given by co-ordinates a,b, and c, ρ is the density of a parcel at a moment t,p is its pressure, ρ 0 and p 0 are the density and pressure of the same parcel under static condition, and x,y,z.Under static condition the parcels obey a static equation ∂p 0 ∂c=−gρ 0 , z=c.One can notice that a nonlinearity of the system of Eqs.(1-5) is associated with the nonlinearity of the inertial terms in the left sides of the Eqs.(1-3) describing parcel's acceleration, the nonlinearity of a mass continuity Eq. ( 4), and the nonlinearity of the adiabatic equation of state (Eq.5).Assume, that there are wave disturbances in the atmosphere that cause small perturbations of parcel's density, ρ 1 =ρ−ρ 0 , or volume V 1 =1 ρ−1 ρ 0 per unit mass, and pressure, p 1 =p−p 0 , relative to their undisturbed values.Denote the components of the displacement vector S(a,b,c,t) of the parcel by x 1 (a,b,c,t), y 1 (a,b,c,t), and z 1 (a,b,c,t), so we can write x=a+x 1 (a,b,c,t), y=b+y 1 (a,b,c,t), z=c+z 1 (a,b,c,t)(6) and where r=(a,b,c), ε αβγ is the anti-symmetric unit tensor of the third rank, δ αβ =1 for α=β and δ αβ =0 for α =β.The sum in Eq. ( 7) is taken over repeating indices α, β and γ .Under incompressible fluid approximation: J =1, which implies a conservation of the density of a fluid parcel along its trajectory.According to the mass continuity equation (Eq.4) the wave disturbances of the volume V 1 are accompanied by the disturbances, J 1 =J −1, of the Jacobian J relative to 1: ∂c , the Eqs.(1-3) can be rewritten in the following form where, according to Eq. ( 5), the volume disturbances are given by The first three nonlinear terms in the expressions for L a , L b and L c describe specifically Lagrangian components of the acceleration, which arise due to transition from Eulerian to the Lagrangian frame.The nonlinearity of the forth (pressure) term in the same expressions is associated with the nonlinearity of the adiabatic state equation (Eq.5).
In case of small relative pressure disturbances, such that 1, the relation ( 13) may be expanded into an infinite power series where c 2 0 = γp 0 ρ 0 is the adiabatic sound speed squared, and via we designated the small terms of order and the terms of higher order.With the use of Eq. ( 14) the Eq. ( 11) takes the following form where is the Brunt-Väisälä frequency squared.
To express pressure terms in Eq. ( 15) via displacements we use Eq. ( 7), which allows us to write the deviation of the Jacobian from 1, J 1 =J − 1, as where J is the sum of the quadratic and cubic (over the displacement derivatives) terms in Eq. ( 7), therefore from Eqs. ( 8) and ( 16) we have ∂x 1 ∂a + Taking the derivative over a of both sides of the Eq. ( 9), and the derivative over b of both sides of the Eq. ( 10), and then adding them with the use of Eq. ( 17) we obtain The Eq. ( 18) allows one to eliminate a pressure term ⊥ ( p 1 ρ 0 ) from the left side of the Eq. ( 15) by applying the operator ⊥ to both sides of this equation.As a result we obtain In what follows we consider internal wave field under incompressible fluid approximation, which requires a formal limit transition c 0 →∞ or V 1 →0 in the Eqs.(8-12), and (14).Such transition implies that the phase speeds of the internal waves considered here are small compared to the sound speed (Gossard and Hooke, 1975).In this case the term ρ 0 V 1 −J = ∂x 1 ∂a+∂y 1 ∂b+∂z 1 ∂c describing linear disturbances of a fluid parcel's volume caused by pressure variations takes the value of the order of p 1 (ρ 0 c 2 0 ) ∼ gz 1 c 2 0 1, therefore the first and last terms in the right side of the Eq. ( 19) can be neglected.As a result Eq. ( 19) takes the form containing only displacement components Assume that the maximum values of the vertical and horizontal displacement components, designated as A c and A a , respectively, are small compared to the corresponding vertical and horizontal length scales, l c and l a , over which a displacement field significantly varies.We consider the ratio µ=A c / l c as a small parameter: µ 1.From a linear continuity equation J 1 ≈ ∂x 1 ∂a + ∂y 1 ∂b + ∂z 1 ∂c ≈0, and from our assumption about the axial symmetry of the displacement field in the horizontal plane we conclude that the following derivatives are small: ∂x 1 ∂a ∼ ∂y 1 ∂b ∼ ∂z 1 ∂c ∼µ 1, therefore the ratio A a l a (or ) is of the order of A c l c .In this case the deviation of the Jacobian (7) from 1, given by Eq. ( 16), contains second-order (∼µ 2 ) and third-order (∼µ 3 ) small terms in the sum J .Now we can estimate the nonlinear terms in the right side of the Eq. ( 20) and compare them with the linear terms in the left side assuming a Boussinesq approximation: l c N 2 g 1, which implies that the mean atmospheric density slowly varies over the vertical length scale l c (Lighthill, 1978).If ω is the characteristic frequency of the internal wave field, 1 and A a ∼ A c l c l a the right side of the Eq. ( 20) takes the value At the same time the left side of the Eq. ( 20) is of the order of ω 2 ( 1 1, therefore the ratio of the terms in the right side of the Eq. ( 20) to the value of the left side ∼ A c l c ∼µ 1.Let us seek a solution of the nonlinear Eq. ( 20) as a series over small perturbations (∼µ ) of the displacement field: z 1 =z 1 +z 1 +..., (similarly, for x 1 and y 1 ), where z 1 =O(µ), z 1 =O(µ 2 ) and etc., to obtain a linear equation for z 1 as a first-order approximation: Then, the second-order equation takes the form where f (z ) is the right side of the Eq. ( 20) expressed through only first-order displacements x 1 =x 1 ,y 1 =y 1 and z 1 =z 1 .Similar equations can be obtained for the higherorder terms in the displacement field.
Our goal now is to consider a process of generation of nonlinear harmonics in a given field of discrete internal gravity waves in the atmosphere and obtain their wave field both in Lagrangian and Eulerian co-ordinate systems.

Second-order Lagrangian displacements induced by internal waves
Below we consider wave-wave interactions based on a wellknown methodology of nonlinear wave theory (see, for instance, Phillips, 1967;McComas and Bretherton, 1977;Craik, 1985;Lighthill, 1978;Rudenko and Soluyan, 1977), but in the Lagrangian frame.Let us present a linear internal wave field in the atmosphere with N=const as a superposition of two plane waves with wave number vectors ) and amplitudes A c,1 , A c,2 .Under Boussinesq approximation: the first-order vertical displacement field z (r,t) (below we omit a subscript 1 for the displacement components) may be presented as where ) 1/2 >0 are wave frequencies, j =1, 2, and k 2 ⊥,j =k 2 a,j +k 2 b,j .Below we consider a two-dimensional case: k b,j =0 and y 1 =0, for which a first-order continuity equation, ∂z ∂c + ∂x ∂a =0, allows us to find a horizontal displacement component Using a linear solution (Eqs.24-25) consider now a secondorder Eq. ( 22) with the initial condition: z (a,b,c,t=0)=0.
Let us seek a solution of Eq. ( 22) in the following form: where A 3 (µt) is a slowly varying amplitude over the time period 2π/ω 3 .Substituting z along with Eq. ( 27) to the Eq. ( 22), and neglecting the terms ∼µ 3 , µ/(k c,3 H ) and 1/(k 3 H ) 2 , we obtain the equation for A 3 (t): then the solution of the Eq. ( 28b) satisfying initial condition The amplitude of the harmonic of difference frequency grows linearly with time where t non is the characteristic interaction time between the two waves, and and periodically varies in time with a period of 4π | ω| .In the latter case |A 3 | reaches a first maximum at a moment t 1 =π/ ω, whereas for the case of resonant interactions it takes a maximum at t=t non in accordance with Eq. ( 29).The ratio of the maximum amplitude under non resonant interactions with | ω|∼ω 3 to that under resonant interactions is the value of order As long as Eq. ( 31) is valid the amplitudes of the nonresonant harmonics are small compared to those of resonant harmonics, therefore wave-wave interactions have a selective character with a nonlinear wave energy exchange taking place mostly between resonant wave modes.However, as a nonlinear parameter µ increases and approaches some finite values smaller than 1, the increasing with µ amplitudes of the non-resonant modes become comparable to those of the resonant modes (according to Eq. 31).In this case the wave energy exchange between wave modes due to non-resonant interactions becomes important as well.Such energy exchange does not have a selective character, since any wave may interact with any other wave to generate a non-resonant harmonic, whose frequency and wave number are not connected via dispersion relation (since ω =0).The importance of the non-resonant interactions between internal gravity waves in shaping of their energy spectrum was pointed out by Phillips (1967), who also noticed that such interactions resemble cascade-like strong interactions between turbulent motions of different scales in stably stratified fluid.Consider first the case of resonant wave-wave interactions.Let designate the aspect ratios of the vertical to the horizontal wave number components of the interacting wave triad as X j ≡ k c,j k a,j , and assume that X 2 j 1 (j =1,2,3).In this case we can neglect by the terms ∼O(1/X 2 j ) in the dispersion equations for each of the three waves, so that the system of the equations describing resonance conditions takes the form (McComas and Bretherton, 1977): After introducing new variables: k a,3 , and eliminating k a,2 k a,3 from the second equation of Eq. ( 32) we obtain the relationship between η and ξ If we choose now a particular set of the vertical wave numbers satisfying the last equation in Eq. ( 32): k c,1 =0.005 rad/m, k c,2 =−0.003 rad/m, k c,3 =0.008 rad/m, then from Eqs. (32-33) obtain: Hence, we can choose the values X 1 =6,X 2 =22,X 3 =8.25 to obtain a set of the horizontal wave number components, k a,1 ≈0.000833rad/m, k a,2 ≈−0.000136 rad/m, k a,3 ≈0.000969 rad/m, and frequencies, ω 1 ≈N/|X 1 |=0.0033rad/s, ω 2 ≈N/|X 2 |=0.0009 rad/s, satisfying resonance conditions (32).Note that a chosen set of the wave numbers for a resonance triad satisfies a Boussinesq approximation, since for all the three waves 1 k c,j H ∼10 −2 , where H ∼8 km.A second-order field for a horizontal component of the displacement,x , can be found from a second-order masscontinuity equation (J 1 =0) for incompressible fluid: where k a,1 , z , x are given by Eq. ( 24) and ( 25), respectively, and the harmonic z has the amplitude A c,3 given by Eq. ( 29).From Eq. ( 34) we obtain Note, that the amplitude of the first term in Eq. ( 35) grows linearly with time, whereas second term remains limited in its absolute value.The expression (29) for the amplitude A c,3 of the difference frequency harmonic was obtained by a small perturbation method, which is valid for the time moments t<t non , for which we can neglect by a reverse effect of the resonant harmonic on the field of the first-order waves.The reverse effect may be taken into account by the method of slowly varying amplitudes (see, for instance, McComas and Bretherton, 1977;Phillips,1967;Craik, 1985), which assumes a slow time variation of the wave amplitudes A c,1 (µt),A c,2 (µt) and A c,3 (µt) over their periods.Using this method one can obtain from Eq. (20) a system of three nonlinear equations for the amplitudes of the wave triad, which shows a conservation of their total wave energy :

Estimates of the amplitudes of nonlinear harmonics in the stratosphere
Let estimate at some moment t<t non the amplitude of the harmonic (31) for the altitudes 20-25 km of the lower stratosphere.At these altitudes the typical rms values of the horizontal velocity fluctuations σ are in the range (0.5-5 m/s) (Fritts and Alexander, 2003), so for N=0.02 rad/s the corresponding rms values of the vertical displacements, ν V =σ/N, are between 25 m and 250 m.If we choose now for the first-order waves the amplitudes A c,1 =65.34 m, A c,2 =59.40 m at some reference height c r =20 km, then from Eqs. ( 27) and ( 29) the amplitude of the harmonic with ω 3 =ω 1 −ω 2 =0.0024 rad/s at a moment t 0 =2.32×10 3 s, equals A c,3 =34.2 m (in this case the coefficient β 1,2 ≈0.41 and the interaction time t non =4.2×10 3 s).Note, that the combinative harmonic with the wave number k 4 =k 1 +k 2 has the corresponding frequency ω 4 =N k a,4 k c,4 =0.007 rad/s that differs from the sum ω 1 +ω 2 =0.0042 rad/s, so this harmonic is a non-resonant one.Its amplitude, designated here as A c,4 , can be estimated from Eq. ( 30) by replacing ω 3 , k 3 on ω 4 ,k 4 , and by taking into account that ω=ω 1 +ω 2 −ω 4 .The coefficient f 1 can be obtained from Eq. ( 21) by changing the signs of k a,2 and k c,2 .In this case the amplitudeA c,4 =27 m at a moment t 0 =2.5×10 3 s, and is comparable with the amplitude of the resonant harmonic.
The measure of nonlinearity of the obtained wave field can be characterized by parameter which is the rms value of the vertical gradient(<(∂z 1 ∂c) 2 >) 1/2 , therefore M 0 is of the same order of smallness as the parameter µ∼ k c,1 A c,1 1 introduced in Sect. 2.
For a chosen set of wave numbers and amplitudes of the wave triad M 0 ≈0.33, and the rms value of the vertical displacements is In the resonance case the first term in Eq. ( 35) is on one order higher in amplitude than the second one, therefore only the first term will be taken into account in the second-order field of the horizontal displacements x 1 =x +x +O(µ 3 ), where x and x are given by Eqs. ( 25) and ( 35), respectively, and via O(µ 3 ) are designated all the nonlinear harmonics of third-order and of higher order.The relative contribution from these harmonics into the total displacement field increases with increasing parameter M 0 .

Vertical profile and spectrum of the wave field in the
Eulerian frame of variables

Resonant interactions
To understand how the obtained Lagrangian wave field is viewed from the Eulerian frame we perform a transformation of variables: x=r+S(r,t), which allows one to calculate the displacement field at some fixed point of space x=(x,y,z) at a moment t.Let designate this field by S E (x,t), and the components of the vector S(r,t) by S a (r,t),S b (r,t), and S c (r,t) instead of x 1 , y 1 and z 1 .Since the displacement S E (x,t) at a point x is caused by a fluid parcel with a co-ordinate r in the unperturbed atmosphere, then S E (x,t)=S(r,t).The latter relationship allows one to express a Lagrangian co-ordinate as r=x−S E (x,t), and obtain the equation for the Eulerian displacement vector S E (x,t): For a given field S(r,t)the solution of the Eq. ( 38) allows one to find the instant vertical profile of the Eulerian displacement vector S E (x=x 0 ,y=y 0 ,z,t 0 ) at a moment t 0 for a given vertical line with the horizontal co-ordinates x=x 0 ,y=y 0 .
In the case of the resonant wave triad, analyzed in Sect.3, the Lagrangian displacements S c (r,t)=z +z +O(µ 3 ) and S a (r,t)=x +x +O(µ 3 ) are given by Eqs. ( 24), ( 25), ( 28) and ( 35).This field will be analyzed in the atmospheric layer of thickness 2h 2H , within which the wave amplitudes A c,j ,(j =1,2,3,...) are considered as constant values taken at some reference height c r of the layer.The substitution of the expressions ( 24), ( 25), ( 28) and ( 35) in the Eq. ( 38) leads to the following system of equations for the vertical and horizontal components, S E,z (x 0 ,z,t 0 ), S E,x (x 0 ,z,t 0 ), of the Eulerian displacement vector S E (x=x 0 ,y=y 0 ,z,t 0 ) at a moment t 0 : where φ 1 , φ 2 , and φ 3 =φ 1 −φ 2 are the initial phases of the interacting waves.
Introducing new non-dimensional variables and performing the following transform of variables: one can rewrite Eqs.(39-40) in the non-dimensional form: For a chosen horizontal co-ordinate x 0 (or R) and a given moment t 0 the first equation in Eq. ( 44) was solved numerically with respect to x .The solution x =x (θ ) was substituted to the Eq. ( 43) to find z =z (θ ), and then Z(θ )=z (θ )−θ from Eq. ( 44).The obtained in parametric form solution of the Eqs.(43-44) allows one to calculate the vertical profiles of the Eulerian displacement components at any time moment t 0 .One of the profiles S E,z (x 0 ,z,t 0 )≡E(z) of the field of resonance wave triad taken at a moment t 0 =2318.4s and at x 0 =3770 m is shown in Fig. 1.Along with the Eulerian field we plotted in Fig. 1. the vertical profile of the Lagrangian vertical displacements S c (a=x 0 ,c,t 0 )≡L(c) of the fluid parcels with the same horizontal co-ordinate a=x 0 in the unperturbed atmosphere and with the different vertical co-ordinates c−c r (in this case L(c) is given by the right side of the Eq. ( 39), where we set S E,x (x 0 ,z,t 0 )=0 and z−S E,z (x 0 ,z,t 0 )=c).
Despite a high nonlinearity of the wave field (M 0 =0.33) the distortion of the Eulerian profile E(z) relative to the Lagrangian profile L(c) is not significant.This fact shows that under resonance conditions the total wave energy is "kept" within a wave triad regardless of the chosen co-ordinate frame.However, the distortion of E(z) becomes significant in the opposite case of non-resonant wave-wave interactions which will be considered below.

Non-resonant interactions
Consider the case, when the interaction between the two Lagrangian waves leads to the generation of the harmonics with combinative frequencies ω 3 =ω 1 −ω 2 , ω 4 =ω 1 +ω 2 and corresponding (related via dispersion equations) wave numbers k 3 and k 4 , for which the synchronism conditions are violated: k 1 −k 2 =k 3 , k 1 +k 2 =k 4 .In this case, a sum in the Eqs.(39-40) should be taken over four terms (j =1,2,3,4).Beside this we have to keep second term in the Eq. ( 35) for the harmonic of the horizontal displacements with ω 3 =ω 1 −ω 2 , and the same term for the harmonic with ω 4 =ω 1 +ω 2 , because these terms can be comparable in amplitudes with the first term in Eq. ( 35).Thus, in the non-resonance case the system of the equations for the Eulerian displacements becomes similar to Eqs. (43-44), but with j =1,2,3,4, and with the following additional term in the Eq. ( 44): where the upper and lower signs correspond to the subscripts 5 and 6, respectively.The chosen set of the wave parameters is typical for the stratospheric altitudes (Eckermann, 1999;Fritts and Alexander, 2003).In this case: B 5 =−5.7,B 6 =−4, r 0 =8.9×10 −5 , ν V =98 m, and the parameter of nonlinearity, M 0 =0.44, shows a strong nonlinearity of the wave field.The calculated instant profiles L(c) and E(z) are shown in Fig. 2.
One can notice a non-sinusoidal " steepening" of the wave crests and troughs in the vertical variations of the Eulerian field E(z), shown in Fig. 2, and a growth of its local gradients as compared to those in the Lagrangian field L(c).Such distortion of the wave field E(z) relative to L(c) is caused by a nonlinearity of the process of advection of fluid parcels induced by the wave field (Chunchuzov, 1996;Eckermann, 1999;CO2).This process is accompanied by a growth of the amplitudes of the high vertical wave number harmonics in the Fourier spectrum of the Eulerian field E(z)relative to those in the Lagrangian frame.At the same time the Lagrangian wave field becomes nonlinear itself while increas- ing M 0 from 0 to M 0 =0.44, because the contribution into this field O(µ 3 ), coming from different nonlinear harmonics, becomes significant and forms a continuous Lagrangian spectral tail due to overlapping of these harmonics.The advection in the Eulerian frame causes additional distortion of the wave field relative to the Lagrangian frame by generating Eulerian spectral tail of higher amplitude than that in the Lagrangian frame.
Thus, the non-sinusoidal distortion of the Eulerian wave field is caused by both the nonlinearity of the wave field in the Lagarangian frame (which distorts the initial Lagrangian field), and the nonlinearity of the transformation from Lagrangian to the Eulerian variables.In the Eulerian frame these nonlinearities cause the non-resonant interactions between original waves and their numerous harmonics, which are open and readily transfer wave energy from the original waves to the high and low vertical wave number harmonics.The non-resonant case as found here significantly differs from the case of resonance wave triad considered above, since the resonance tends to keep energy within a wave triad and prevents its distortion.

Vertical wave number spectral tail generated by a nonlinearity of the wave field
The nonlinear distortion of the Eulerian vertical profile E(z) is accompanied by enhancement of the amplitudes of high vertical wave number harmonics of given two Lagrangian waves in the Fourier spectrum of E(z).These harmonics are generated due to numerous non-resonant wave-wave interactions, and form a high vertical wave number tail in the power spectral density (PSD) of E(z) defined as (Bendat and Piersol, 1967): where F (k z ) is a Fourier transform of E(z) calculated below by using MATLAB toolbox, and l is the length of the realization of E(z).While increasing the parameter of nonlinearity M 0 from M 0 =0.05 (Fig. 3a) to M 0 =0.44 (Fig. 3b) the amplitude of the tail in the Eulerian spectrum increases as well.For small wave amplitudes, M 0 =0.05, there is a very weak difference between Lagrangian and Eulerian spectra, although the appearance in the Eulerian spectrum of the peak corresponding to the second-order harmonic with k c,3 =k c,1 −k c,2 =0.008 is already noticeable (see continuous line in Fig. 3a).The amplitudes of the second-order harmonic and of higher order harmonics in the Lagrangian and Eulerian spectra (shown in Fig. 3b by the dotted and continuous lines, respectively) gradually increase with increasing M 0 .However, the most rapid increase of these harmonics takes place in the Eulerian spectrum, whose high-wave number tail approaches at M 0 =0.44 the straight line corresponding to the observed power law decay βk −3 z , β=0.1.While changing the location x 0 of the profile E(z), the slope of the spectral tail also changes due to horizontal nonhomogeneity of the wave field, but an average power law decay of the tail remains close to the βk −3 z -power law (Chunchuzov, 2008).It is important to note that our analysis of the internal wave field was based on the approximate solutions of the motion equations for nonviscous fluid (Eqs.1-5).These solutions do not take into account the existing mechanisms of wave energy dissipation in the realistic atmosphere (Gossard and Hooke, 1975) that limit an infinite growth of the local wave field gradients with increasing value of M 0 .In the absence of dissipation a strong nonlinear steepening of the wave crests in the Eulerian profile E(z) for M 0 =0.44 (Fig. 2) leads to the extremely high values of the vertical gradients ∂S E,z /∂z , which can exceed 1 at certain altitudes, as seen in Fig. 4a.Beside this, the calculated for the secondorder displacement field vertical variations of the Jacobian J 1 =J −1 (see Fig. 4b), have the rms values of about 0.2, and at certain altitudes may reach maximum values of about 0.6.This means that the Lagrangian displacement field containing second-order terms (∼M 2 0 ), obeys the approximation (34) of fluid incompressibility, J 1 =J −1=0, only with the accuracy of the terms of order M 2 0 .The third-order terms (∼M 3 0 ) like ∂z ∂c ∂x ∂a − ∂z ∂a ∂x ∂c +..., and higher-order terms contribute to the deviation J −1 from incompressibility, which increases with increasing value of M 0 .Therefore, the accuracy of the incompressible fluid approximation, which was initially assumed in the paper for small enough values of M 0 , decreases with an increase of M 0 .
The high-order terms in the displacement field contain the harmonics with the combinative frequencies pω 1 ±qω 2 (p and q are integer numbers) greater than the buoyancy frequency N. As a result the temporal variations of the Eulerian displacement field become also nonsinusoidal and contain a nonlinear high-frequency tail (∼ω −3 ) in their frequency spectrum (Chunchuzov et al., 2006).Thus, a highly nonlinear displacement field may vary in time so fast that during short characteristic time scales τ <2π/N the sound travels on the distances c 0 τ that are less than the vertical length scales 2π/k z of the field's spatial variations.For such short time intervals τ the pressure of adiabatically displaced fluid parcel is unable of adjusting to the hydrostatic atmospheric pressure, since this parcel "meets" a background, which is already disturbed by fluid compressibility.Therefore, there is no adiabatic constrain imposed by KS06 on the volume variations J −1, because this constrain was based on the assumption about the existence of the parcel's pressure adjustment to the hydrostatic background atmospheric pressure.Such assumption is valid for only slow displacements with small amplitudes, but not for nonlinear displacements with high amplitudes and fast temporal variations.
To prevent the extremely large deformations of fluid parcel's volume J 1 =J −1, arising in the ideal fluid under strong nonlinearity of the wave field, we have to limit the spatial gradients of the obtained Eulerian and Lagrangian displacement fields by taking into account the dissipation of wave energy.Otherwise, the further increase of the parameter M 0 from 0.44 to 0.52 leads to the local absolute values of J 1 exceeding 1 at certain altitudes, what means the nonphysical 100% relative changes of the fluid parcel's volume and the invalidity of the Lagrangian to Eulerian variable transformation.In the latter case the vertical profile E(z) becomes discontinuous and multivalued in the thin local regions of space (shown in Fig. 5a by arrows) due to unlimited nonlinear distortion of the wave field.Such multivalued Eulerian profiles were earlier found by KS06 and Klaassen (2009) in their kinematic advection models for a non-dissipative wave field comprised of a superposition of seven or more linear Lagrangian gravity waves.These models extended earlier Eckermann's model (Eckermann, 1999) by taking into account compressible gravity waves.
It was found in KS06 that for the typical wave amplitudes in the middle atmosphere required for producing a broad Eulerian tail with a slope of −3, the Lagrangian to Eulerian transformations are overwhelmingly singular, with multiple parcels frequently occupying the same physical points in space.These singularities are induced by stretching deformation fields which form during the superposition of sinusoidal waves with nonparallel wave vectors.Such deformation fields were assumed in K09 to be unstable with respect to three-dimensional vortices, therefore ". . . the saturated middle atmosphere wave fields are frequently accompanied by small-scale turbulent eddies".Thus, the instability proposed in K09 leads to the wave breaking processes which are seen (Fig. 5a) to occur in the thin (as compared to the vertical wave lengths of the main Lagrangian waves) atmospheric layers with the high negative vertical gradients of the vertical displacements, ∂S E,z /∂z<−1.These layers have high positive vertical potential temperature gradients ∂θ/∂z=∂θ 0 /∂z(1−∂S E,z /∂z), where θ 0 is unperturbed potential temperature.
Beside the thin layers with singularities there are another layers, where the vertical gradients of the vertical displacements are positive and reach their critical values, ∂S E,z /∂z≈1, at which wave induced convective instability switches on and generates turbulence (Gossard and Hooke, 1975).In these layers the wave-perturbed potential temperature gradient ∂θ/∂z becomes close to zero, but each wellmixed turbulent layer is confined between two thin layers with high positive potential temperature gradients or absolute temperature inversions (Whiteway et al., 1995).The convective instability, however, arises for higher values of M 0 than the condition ∂S E,z /∂z<−1.
Thus, wave breaking processes of different types generate turbulent eddies of different scales confined within the local regions of space, where ∂S E,z /∂z reaches the values of order 1.These regions are the main sinks of wave energy, whose dissipation stabilizes the nonlinear growth of the wave field gradients by some finite values.In CO2 the vertical length scale l c =2π/m c of these regions was estimated for a large number (n 1) of random waves by limiting the rms values of the gradient ∂S E,z /∂z by 1, where m c is the upper vertical wave number of the spectral tail.For the stratosphere the estimated value of l c was in the range (10-100) m, which was close to the characteristic scale found in the spectra of star scintillations measured from space stations (Gurvich andChunchuzov, 2003, 2005).This scale separates the two regions of star scintillation spectrum associated with the anisotropic and isotropic temperature fluctuations in the stratosphere.Taking l c ∼50−100 m, we have limited the calculated spectrum in Fig. 3b by the corresponding critical vertical wave number m c .
There are two competitive processes that are shaping the Eulerian spectral tail.From one side, a strong nonlinearity of the wave field leads to the non-resonant wave-wave interactions that become significant for the vertical wave numbers k z >ν −1 V .The characteristic vertical wave number m * =1/(2 1/2 ν V )=0.007 rad/m above which the non-resonant interactions form spectral tail is indicated in Fig. 3b (the corresponding vertical scale L * =2π/m * ≈900 m).These interactions transfer wave energy from the wave numbers of the most energy contained waves, excited by wave sources, toward higher vertical wave numbers k z >m * , although some energy transfers to the lower wave numbers as well (as seen in Fig. 3b).From the other side, at high vertical wave numbers k z approaching a critical wave number m c m * , the dissipation of the wave energy due to wave breaking processes becomes important.These processes generate the turbulent eddies of different vertical scales 2π/k z <l c within the thin layers co-existing with the stably stratified layers of larger vertical scales, The arising turbulent diffusion smoothes high gradients of the wave field (Gossard and Hooke, 1975;Weinstock, 1985;Whiteway et al., 1995), thereby preventing extremely large volume deformations of the fluid parcels.At the same time the arising turbulent layers become confined by the thin stable layers with absolute temperature inversions (Whiteway et al., 1995).The resemblance of the stable layers of high positive potential temperature gradients with the narrow "sheets" observed in high resolution temperature profiles (Dalaudier, 1994) was first noted by Eckermann (1999).
5 Eulerian vertical wave number spectrum of ducted gravity wave mode.

Nonlinear distortion of ducted wave field
Although, the most part of the gravity wave spectrum in the atmosphere is composed of freely propagating waves (unlike the ocean), there exist stably stratified layers in the atmosphere where the trapped waves constitute a high-frequency part of the spectrum with periods usually less than 10-15 min (see, for instance, Böhme et al., 2004).To show the generality of the proposed shaping mechanism we will find below the Eulerian vertical spectrum of the displacement field induced by only one standing gravity mode in the atmospheric layer of 2h thickness.Within this layer we take the buoyancy frequency N=const for −h≤c≤h and N 1 N outside this layer.For ducted wave the Lagrangian displacements may be written as S a (r,t) = XAcos(mc)sin(ωt where 2π m is the vertical period of standing wave, A is the modal amplitude, ω=N |k a | (k 2 a +m 2 ) 1/2 and ϕ are the modal frequency and phase, respectively, X=m/k a , and via O(µ 2 ) are denoted the second-, third-and higher-order fields generated in a given Lagrangian field of ducted mode.For ducted wave (Eqs.47-48) the system of the Eqs.(38) for the Eulerian displacement field S E =(S E,x ,S E,z ) taken at x=x 0 and at a moment t 0 becomes where The Eqs. (49-50) may be solved similarly to those derived in Sect. 3 for the wave triad (Eqs.43-44).For the time moments t 0 and horizontal co-ordinates R 0 satisfying a condition: ωt 0 +k a R 0 +φ=2πn,n=0,±1,±2,..., the horizontal displacement component x =0, therefore x 0 =a.In this particular case Eq. ( 49) takes a simple form which is similar to the solution of the equation describing acoustic wave of finite amplitude (Rudenko and Soluyan, 1977).The solution of Eq. ( 52) may be written in the form It is seen from Eq. ( 53) that for small relative amplitudes of gravity mode (M 1) the vertical profile z =z (Z,t 0 ,R 0 ) of the Eulerian vertical displacements is almost of the same sinusoidal form as that for the Lagrangian vertical displacements with the horizontal co-ordinate a=x 0 (where S E,x =S a =0).While increasing M up to 0.9 the Eulerian profile S E,z (z,t 0 )=Az (shown by continuous line in Fig. 6a) is distorted due to advective nonlinearity so that the points on the sinusoidal profile S c (c,t 0 ) (dotted line in Fig. 6a) of the Lagrangian standing wave (Eq.52) are differently displaced with respect to the nodes of the standing wave.Such displacements grow with the local absolute values of S c (c,t 0 ), and this leads to the non sinusoidal "steepening" of the Eulerian profile S E,z (z,t 0 ) relative to the profile S c (c,t 0 ).One period of the standing wave (Eq.52) is shown in Fig. 6b.This fragment illustrates the vertical displacements S c (c,t 0 ) of the selected four equidistant fluid parcels with Lagrangian co-ordinates c 1 ,c 2 ,c 3 ,c 4 .The parcels are instantly displaced to the points with the vertical coordinates x 1 ,x 2 ,x 3 ,x 4 so that x i =c i +S c (c i ,t 0 ),i=1,..,4, and S E,z (z=x i ,t 0 )=S c (c i ,t 0 ).Such displacements condense the positions x 2 ,x 3 ,x 4 of the parcels 2, 3, 4 as compared to their undisturbed positions c 2 ,c 3 ,c 4 , and, at the same time, increase the distances between parcels 1 and 2. The non sinusoidal distortion of the Eulerian profile S E,z (z,t 0 ) is accompanied by an increase of the local vertical gradients ∂S E,z /∂z, and by generating of the high vertical wave number harmonics of the main wave (Eqs.47-48) in the Fourier spectrum of S E,z (z,t 0 ).The overlapping spectral picks corresponding to the different harmonics form a nonlinear tail in the Fourier spectrum.

Vertical wave number spectrum of ducted wave field
The vertical wave number spectra for the displacement fields S E,z (z,x 0 ,t 0 ) at x 0 =0 and ϕ=π/2 are shown in Fig. 7 for different time moments t 0 .These spectra are calculated for M=0.9, for which the parameter of nonlinearity M 0 =M/2=0.45.One can find well-pronounced spectral peaks corresponding to the second-, third-and higher-order harmonics of the main mode with m=0.0075 rad/m.Despite a "brush-like" form of the obtained spectra their power law decays are close to the observed βk −3 z -form, and such agreement supports our assumption about a nonlinear origin of the observed spectral tail.
A further growth of M up to 1 (or M 0 up to 0.5) leads to the multivalues in the profile S E,z (z,x 0 ,t 0 ) near the points Z p =π+2πp, (p-integer), x =0 (or S E,z =0), and z =0(or S E,x =0).
At these stagnation points the gradient of the displacement field (Eq.52) In the realistic atmosphere as known the turbulence generated by wave breaking processes smoothes high gradients of the wave field due to turbulent viscosity (Gossard and Hooke, 1975), thereby preventing the arising of discontinuities and of extremely large fluid parcel deformations.However, such smoothing as mentioned above occurs only in the local regions of the wave field which have relatively small vertical sizes l c ∼(10-100) m, whereas wave field vertical variations of larger scales, l c <2π/k z <2π/m * , are caused by various nonlinear harmonics generated due to non resonant wavewave interactions.These interactions lead to the cascade-like wave energy transfer over wave numbers of the spectral tail towards high vertical wave numbers k z ∼m c , at which wave energy transfers to the energy of turbulent eddies.
The role of nonlinearity and dissipation in shaping of the vertical profile of nonlinear gravity wave field has certain similarity with that in shaping of so-called N-wave in the nonlinear acoustics (Rudenko and Soluyan, 1977).The latter is the exact solution of the Lagrangian wave motion equations for ideal fluid (Zarembo and Krasilnikov, 1965), which describes a nonlinear distortion of initially sinusoidal wave profile with increasing distance from a plane acoustic source.The points on the acoustic wave profile propagate with the different local speeds proportional to the fluid parcel velocities in the acoustic wave.This causes a steepening of the wave form with increasing distance from a source so that wave profile in the nonviscous fluid ultimately becomes discontinuous and multivalued near wave front.The "curing" of the profile from multivalues can be achieved by taking into account a molecular viscosity and thermal conductivity.At large Reynolds number, Re 1, the competition between the nonlinear distortion of the wave profile and dissipative smoothing of high gradients at the wave front establishes Nlike wave form at some distance from a source.It is important to note that dissipation takes place mostly within a thin wave front with a thickness ∼1/Re, whereas the rest part of the N-wave remains distorted due to nonlinearity of the acoustic wave field.For highly nonlinear gravity waves considered here a smoothing of the large wave field gradients in the middle atmosphere is caused mostly by turbulent viscosity.

Analytic solution for a periodic ducted wave field
Consider now a limiting case, when the thickness 2h of the atmospheric layer tends formally to infinity.For −∞<z<∞ a vertical displacement field x =F (Z) obeying Eq. ( 52) is a periodic process with a 2π/m period.In this case we can find a solution of the Eq. ( 52) and its discrete power spectrum S(k z ) in analytic forms (Appendix A).Indeed, a solution of the Eq. ( 52) may be presented as a Fourier series: S E,z (z)=Az =A ∞ n=1 b n sinnmz, with the coefficients b n (M)= −2(−1) n J n (nM) nM for n=1,2,. . ., (given by Eqs.A1 and A6).Using a correlation function of the periodic displacement field S E,z (z) (formula A9) and the corresponding discrete power spectrum defined by Eq. (A10) one obtains the following expression for one-sided power vertical wave number spectrum S(k z ): It is important to note that the intensities c 2 n of the spectrum (Eq.54) at discrete vertical wave numbers k z =nm depend only on the non-dimensional amplitude M=mA (or M 0 ) and the harmonic's number n.For n=0 the spectrum has a zero intensity, since the mean (over period) value <S E,z (z)>=0.For growing values of M a spectral intensity c 2 n as a function of k z =nm is shown in Fig. 8.While increasing M from 0.7 to 0.9 the discrete intensities of the spectral tail with n>1 increase up to the values, whose power law decay with increasing n becomes close to n −3 (a straight line 0.4/(m 2 n 3 ) in Fig. 8) within some range of n.However, at high n they show a more rapid decrease with increasing n than n −3 .The vertical wave number range of the spectral tail is limited by critical wave number The 2-D case considered above demonstrates the simplest examples with only two propagating waves and one ducted wave, when the nonlinear mechanism works.The applicability of the 3-D case in the real atmosphere is shown below for an ensemble of gravity waves with arbitrary amplitudes and phases.
6 Eulerian spectrum of a random ensemble of internal gravity waves in the atmosphere

Relationship between Eulerian spectra and Lagrangian correlation functions of the random waveinduced displacements
For the case of high number n of waves (n 1) with randomly independent amplitudes and phases, previously analyzed by AJ89, Chunchuzov (1996), Chunchuzov (2001)  equilibrium internal wave spectrum is the statistical approach based on the transformation (6) from Lagrangian to Eulerian variables: x=r+S(r,t), S E,j (x,t)=S j (r,t),(j =1,2,3), which can be also derived through the delta function δ(x) as follows: S E,j (x,t)= d 3 rJ S j (r,t)δ(x−y), y=r+S(r,t), r=(a,b,c),x=(x,y,z), and J is the determinant of the Jacobian of transformation defined by Eq. ( 7).Using a Fourier transform δ(x−y)=(2π) −3 d 3 kexp[ik(x−y)] one obtains the relationship between Lagrangian and Eulerian vertical displacements (AJ89, C96, CO2, H01): Based on the assumption about the random displacement field S(r,t) as being a statistically stationary, homogeneous and Gaussian field the transformation (55) allows one to obtain a relationship between the 3-D Eulerian wave number spectrum S E (k) and the correlation functions of the Lagrangian displacement field (CO2): In general case of J =1 the linear, quadratic and cubic (over the displacement gradients) terms in Eq. ( 7) describe the deformations of a fluid parcel due to its compressibility.They lead to the additional terms in G s (R,k) like those produced by the cubic terms: At high wave number values only small-scale variations of the displacement field S(R,t)significantly contribute to the spectrum (Eq.56).If at small scales R i a structure function 2D ij (R) is of quadratic form: then at high k the 3-D spectrum (Eq.56) has a k −5 power law decay, and the corresponding vertical wave number spectrum takes a universal k −3 z -form at high k z .This general result was obtained in CO2 for any Lagrangian displacement field S(R,t) having a quadratic structure function (57).Such structure of the field S(R,t) is typical for both linear and nonlinear gravity waves.However, for estimating the coefficients C ij mn in Eq. ( 57) the field S(R,t) was chosen in CO2 as a superposition of weakly interacting gravity wave modes.

3-D and 1-D spatial spectra of the wave-induced displacements
The nonlinearity of the Lagrangian motion equations (1-5) causes the interactions between gravity waves so that some equilibrium wave energy distribution among Lagrangian wave modes forms due to the energetic balance between the nonlinear energy transfer from the characteristic horizontal −1 0 ) and vertical (m −1 0 ) scales of the source spectrum toward smaller scales, and the wave energy dissipation at small vertical length scales.It was shown in AJ89, HO1, and CO2 that the form of the Lagrangian wave energy distribution E(k ⊥ ,m) between modes does not affect the asymptotic form of the 3-D Eulerian spectral tail at high vertical wave numbers, therefore E(k ⊥ ,m) was chosen by CO2 in the form: For 2 1/2 |k z |ν V 1, (ν V is the rms value of the vertical displacement field), and for the low horizontal wave numbers 1, the 3-D spectrum (Eq.56) takes the asymptotic form: where χ=m 0 /k 0 ≈ν 1 /ν V , ν 1 is the rms value of the horizontal displacement components, a 0 =−C 3333 /16 and e 0 =−C 3311 /8 are the coefficients proportional to the mean square values of the vertical and horizontal gradients of the vertical displacement field S c (r,t), respectively.Using Eq. ( 58) the quantities a 0 and e 0 may be expressed through the parameters of nonlinearity M≡m 0 ν V and χ: a 0 = M 2 8 , e 0 ≈ a 0 χ 2 (the parameter Mwas introduced in CO2 and equals M 0 /2 1/2 ).They characterize the degree of nonlinearity and anisotropy of the wave field, respectively.
The obtained 3-D spectrum (Eqs.59-60) shows the existence of highly anisotropic spatial inhomoigeneities in the displacement and temperature fields, whose anisotropy depends on the value of e 0 1.This spectrum was used by Gurvich andChunchuzov (2003, 2005) for explaining the spectra of stellar scintillations observed from space, and by Ostashev et al. (2005) for modelling of a scattering of the acoustic waves from anisotropic wind speed and temperature inhomogeneities caused by internal wave field in the atmosphere.
After integrating Eq. ( 59) over 2π k ⊥ dk ⊥ from 0 to ∞ we obtain a 1-D vertical wave number spectrum of the vertical displacements S 1E (k z ): where β is given by Eq. ( 60), m * = 1 2 1/2 ν V = N 2 1/2 σ is the characteristic vertical wave number, above which the nonresonant wave-wave interactions become important and form a universal spectral tail (Eq.61), σ is the rms value of the Lagrangian velocity fluctuations, and m c =m * exp( 1 β ) is the critical vertical wave number, at which wave energy dissipates due to wave instabilities.
As seen from Eq. ( 60) the amplitude β of the spectral tail (Eq.61) increases rapidly with increasing a 0 (or M), but starting from a 0 ≈0.01 such increase significantly slows down (Fig. 9a) so that β reaches a broad maximum of about 0.22 at a 0 ≈0.012 (M∼0.315).Therefore, a further increase of M does not change the amplitude of the tail.Moreover, the growth of M is limited itself by the wave-induced instabilities that prevent the parameter M to reach the value M≈0.5 for which <(J −1) 2 >=1.In the nonviscous fluid the latter condition leads to the discontinuities in the wave field profile and the divergence of the interaction potential energy of the internal wave field as shown in AJ89 and CO2 (p. 1771 ).
Note that the amplitude of the spectral tail β reaches a maximum of about 0.2 and becomes saturated for the values M=0.3-0.4,for which the variance of the Jacobian, <(J −1) 2 >, being the value of order M 4 (CO2, p. 1771), remains less than 1.In the local regions of space, where a highly nonlinear wave field breaks into turbulence, the turbulence viscosity is expected to smooth extremely high gradients of the wave field.Such smoothing also decreases the value of <(J −1) 2 > as compared to that in the nonviscous fluid.Thus, the obtained k −3 z spectral tail is saturated as a result of a combined effect of strong nonlinearity of the wave field within wave number range (Eq.61) and dissipation of the wave energy at high values of k z near the critical wave number k z ≈m c .The same amplitude dependence with growing M 0 we obtained earlier for the spectrum of a few discrete waves (shown in Fig. 5).Hence, a process of nonlinear saturation of the amplitude of spectral tail takes place for both a few discrete waves and for a high number of random waves with a broadband wave number spectrum.
At low wave numbers the 3-D spectrum SE (k ⊥ ,k z ) significantly depends on the power spectrum of the random internal wave sources, and is weakly affected by the nonlinear effects.Therefore, for 2(k 2 ⊥ ν 2 1 +k z ν 2 V ) 1 the forms of the 3-D Lagrangian and Eulerian displacement spectra are only slightly differ from each other (see AJ89 and CO2) .In this case the 3-D spectrum SE (k ⊥ ,k z ) can be chosen in the form given by Eq. ( 58), where C is some constant.For k 2 z m * 2 the 1-D vertical wave number spectrum of the vertical displacements, S 1E (k z ), can be obtained by integrating Eq. ( 58) over 2π k ⊥ dk ⊥ from 0 to ∞.
The low wave number part of the spectrum S 1E (k z ) is connected at the intermediate wave number k z ∼m * with a high vertical wave number tail (Eq.61) by some transitional portion of the spectrum, whose exact shape is unknown to us.To find an approximate shape of this spectral portion we take into account the continuity of the spectrum at k z =m * .Equating at k z =m * the low-wave number part of the spectrum S 1E (k z ) and its high-wave number tail we can find the constant C. Taking also into account that where the ratio m * /m 0 depends on β.For M=0.34 the spectral tail is under saturation condition, for which β=0.22 and m * /m 0 ≈2.For this case the spectrum (Eq.62-63) is shown in Fig. 9b.for the two rms vertical displacements, ν V =244 m and ν V =623 m, corresponding to the two different altitudes 15 km apart.The increase of ν V with altitude due to decrease of the mean atmospheric density leads to the decrease of the wave number m 0 of the spectral maximum along with the characteristic wave numbers m * and m c .Such increase, however, does not change the amplitude of the tail (Eq.63), whose characteristic (outer) vertical scale 2π/m * increases from 2167 m to 5712 m along with an increase of the wave breaking scale 2π/m c from 23 m to about 59 m.It is necessary to note that a generation of the k −3 z -tail in the Eulerian frame of variables is not a purely kinematic effect associated with a nonlinear transformation from Lagrangian to the Eulerian frame of variables: the advective nonlinearity generates a saturated tail for only those values of M 0 for which a nonlinearity of the Lagrangian wave field becomes important as well.

Conclusions
In this paper a nonlinear shaping mechanism was studied for the vertical wave number spectrum of the field composed of two propagating gravity waves or one ducted gravity wave.This mechanism is similar to that previously proposed in CO2 for an ensemble of large number of waves with randomly independent amplitudes and phases.In both CO2 and present paper we took into account a strong nonlinearity of the gravity wave field for the wave amplitudes and wave numbers typical for the real atmosphere.This was done by solving the nonlinear fluid motion equations in the Lagrangian frame and by applying a variable transformation to the Eulerian co-ordinates.Such an exact transformation allowed us to take strictly into account the advection of fluid parcels induced by the wave field itself.The different advection of different fluid parcels as shown leads to the "steepening"of the wave crests and troughs when transferring from Lagrangian to the Eulerian frame of variables.Such non sinusoidal distortion of the wave field vertical profile is accompanied by generating of a high vertical wave number spectral tail in its Fourier spectrum.The 3-D form for this tail was used in the paper to obtain a new form (given by Eqs.62-63) for the 1-D vertical wave number spectrum in the broad range of wave numbers.
The amplitude of the spectral tail increases with the parameter of nonlinearity M 0 up to some broad maximum or saturation for which the 1-D spectral tail takes a form close to βk −3 z regardless of whether wave sources excite in the Lagrangian frame only two discrete internal waves, one ducted wave or a random ensemble of waves with a broadband wave number spectrum.This was also confirmed by an analytic solution found for a discrete vertical wave number spectrum of the periodic (over z) standing wave of finite amplitude.A further increase of M 0 , and of the spatial gradients in the Eulerian displacement wave field, is limited by their threshold values, at which the wave breaking processes generate turbulent eddies in certain local regions of space.
Thus, on reaching M 0 certain values (about 0.3-0.4)there exists some vertical wave number range within which the wave energy transfer toward high vertical wave numbers is balanced by a sink of the wave energy due to wave breaking processes at some critical wave number.Within this range a vertical wave number spectrum of the Eulerian vertical displacements takes a universal βk −3 z -form with the theoretical value β=0.22 that lies within its observed range (0.1-0.3).

Fig. 1 .
Fig. 1.One of the vertical profiles E(z) (continuous line) of the Eulerian vertical displacement field of the resonant wave triad calculated for M 0 =0.33, x 0 =3770 m and t 0 =2318.4s.Also shown are vertical displacements S c (a=x 0 ,c,t 0 )=L(c) of the parcels with the same Lagrangian horizontal co-ordinate a=x 0 and different vertical co-ordinates c−c r (dashed line), c r is some reference co-ordinate c taken inside a given atmospheric layer.

Fig. 3 .
Fig. 3. Evolution of the vertical wave number spectra for the vertical profiles E(z) (continuous line) and L(c) (dotted line) with growing parameter of nonlinearity M 0 .(a) M 0 =0.05,νV =12 m.(b) M 0 =0.44,νV =98 m; the corresponding profiles are shown in Fig. 2. The straight dashed line corresponds to the observed power law decay βk −3 z , β=0.1.The wave numbers m * and m c are indicated by arrows, where m * =1/(2 1/2 ν V ) is the lower characteristic wave number above which the non-resonant interactions form spectral tail, and m c is the critical vertical wave number, at which wave energy dissipates due to wave-induced instabilities.

Fig. 4 .
Fig. 4. Vertical profiles of the vertical gradient of the Eulerian displacements, ∂S E,z /∂z, and of the relative volume deformations of fluid parcels, J 1 =J −1, for the case shown in Fig. 2. (a) ∂S E,z /∂z vs. vertical Eulerian co-ordinate z−c r .Vertical straight line corresponds to the value ∂S E,z /∂z=−1.(b) J 1 =J −1 vs. vertical Lagrangian co-ordinate c−c r .

Fig. 6 .
Fig. 6.Vertical profiles of the vertical displacement field of ducted gravity mode in Lagrangian (dotted line) and Eulerian (continuous line) frames,M=0.9;(a) distortion of the profile S E,z (z,t 0 ); (b) fragment of S E,z (z,t 0 ) illustrating the distortion of the profile due to displacements S c (c,t 0 ) of the selected four equidistant fluid parcels with the Lagrangian co-ordinates c 1 ,c 2 ,c 3 ,c 4 .
m c =mn c , at which the variance of the displacement gradient, <( ∂S E,z (z) ∂z ) 2 >= n c n=1 m 2 n 2 c n 2 = n c n=1 2J 2 n (nM), reaches a value of about 1.The critical vertical wave numbers m c = mn c above which wave field breaks into turbulence are indicated by vertical arrows in Fig. 9.The corresponding minimal vertical scales, 2π/m c , are of about 34 m for M=0.87, and 93 m for M=0.9.It is of importance that the analytically obtained discrete spectrum shows the same evolution of its tail with increasing parameter M as numerically estimated continuous spectral tail of ducted gravity mode in the layer of finite thickness −h<z<h.

Fig. 9 .
Fig. 9. Vertical wave number spectrum of a random internal wave field, given by Eqs.(62-63), and its amplitude β.(a) β as a function of the parameter a 0 ≡M 2 /8 characterizing the degree of nonlinearity of the wave field.(b) Model spectrum S 1E (k z ) calculated for two altitudes, which are 15 km apart.