Weakly Nonlinear Stage of Instability Development in a Sharply Stratified Shear Flow with an Inflection–Free Velocity Profile

The distinctive property of the class of shear flows under study is that in a large part of the instability domain the phase velocities of waves are so close that their individual critical layers merge into a common one. Throughout a weakly nonlinear stage of perturbation development, this is the layer in that the most intensive and diverse wave interactions operate which determine scenario of perturbation evolution. Analysis of these interactions allows us, first, to reveal two stages of the evolution, three-wave one, when three-wave interactions dominate, and post-three-wave when numerous nonlinear interactions of different orders come into play, and, second, to determine which of the higher-order interactions are competitive. On this basis, we have found the structure of nonlinear evolution equations, substantiated that the nonlinear growth of wave amplitudes is explosive, and calculated growth indexes for both nonlinear stages. It is found that during the three-wave stage the most rapidly growing are low-frequency waves whereas at the next stage the growth of high-frequency waves is accelerated, and to the end of the weakly nonlinear stage all the waves have amplitudes of the same order. The results obtained are illustrated by numerical calculations for some ensembles of waves.


Introduction
In astrophysical and geophysical fluid dynamics, high-Reynolds-number shear flows are often sharply stratified, i.e. they are stratified in such a manner that the vertical scale of density variation ℓ is much smaller than that of velocity shear, Λ [1][2][3][4][5]. In these circumstances, a key role in the flow development (in time and/or space) is played, as a rule, by Holmboe's instability [6] which attracts a considerable interest of researchers. Now their efforts are concentrated in two main directions. First, the linear theory of instability is developed [7][8][9][10][11] which describes the properties of Holmboe waves at the initial stage of growth. Secondly, both laboratory experiments and direct numerical simulations are employed to explore the strongly nonlinear stage of their evolution. In this way, a number of interesting and important results were obtained concerning the structure and spectra of Holmboe waves, as well as the mixing and vertical transport due to these waves (see [12] and references therein).
However, it would be appropriate to distinguish yet another stage of unstable disturbances growth which is transitional between the linear and strongly nonlinear ones. It is at this stage that all the structures and features of Holmboe wave field observed in numerical simulations and laboratory experiments appear and begin to develop, and therefore its studying is crucial for understanding the mechanisms of the instability. Usually, such an intermediate stage is described by weakly nonlinear theory, and the aim of this paper is to study a weakly nonlinear stage of Holmboe's instability development. To do this, let us consider a class of flows in which Holmboe waves are the only possible waves. In particular, among these are flows of an ideal incompressible fluid which are similar to a boundary layer with an imbedded thin pycnocline. Their velocity V x = U (z) increases monotonically upwards, from zero at the bottom (z = 0) to some U 0 as z → ∞, whereas the density ρ 0 (z) steadily decreases from ρ 1 to ρ 2 in a transition layer of thickness ℓ ≪ Λ (Figure 1). In what follows, we shall deal with such flows.
It is well known that the key role in the development of an instability in high-Reynolds-number shear flowsduring both the linear and nonlinear stages -is played by the resonant interaction between the wave and the flow [13][14][15] in so called critical layer [16,17] (hereafter referred to as CL), i.e., in a narrow vicinity of the critical level z = z c on which the flow velocity is equal to the wave phase velocity, U (z c ) = c r ≡ Re c. Owing to the resonance, an intensive momentum and energy exchange between the instability wave and fluid particles inside CL takes place, and the actual magnitude of perturbation is there much greater than in any other part of the flow. At the linear stage, perturbation can be regarded as a superposition of independent eigenoscillations, each with its own phase velocity and CL. With the growth of the amplitudes, however, this approximation ceases to be valid, and a weakly nonlinear stage of evolution begins. During this stage, the perturbation can still be represented as an ensemble of waves arranged in almost the same manner as in the linear approximation, but it is impossible to neglect their interactions. And it is the CL where the perturbation magnitude is the greatest and the most intensive wave interactions take place [17,18]. Hence, a significant contribution to the perturbation evolution can be made only by those interactions which involve the waves with a common critical layer, i.e., with nearly the same phase velocities (socalled phase-locked modes, see [19,20], for example).
The linear theory of Holmboe waves in flows of the class under study [21,22] demonstrates that their properties differ essentially from those of waves in homogeneous shear flows. For a brief review of them let us introduce dimensionless variables so that (prime denotes the derivative in z) and write the squared buoyancy frequency as Ω 2 (z) ≡ −gρ ′ 0 /ρ 0 = J n(z) where g is the gravitational acceleration, and J is the bulk Richardson number. The normalized buoyancy profile n(z) is localized in a layer of thickness O(ℓ) centered at z N = O(1), and has a single maximum.
In the Boussinesq approximation, linearized hydrodynamic equations for a single harmonic of planar (i.e., independent of y) disturbance reduce to Taylor-Goldstein equation [16] with boundary conditions where w(z) exp[i(kx − ωt)] is the vertical velocity component, and c = ω/k is (complex-valued) phase velocity.
Since Ω 2 (z) decreases rapidly above the pycnocline, each wave has a reflection point at a certain height. Therefore, the flow can be regarded as a horizontal waveguide, and its eigenoscillations as guided modes. Waves overtaking the flow are neutrally stable (Im c = 0) and characterized by wavelength and by the number m = 0, 1, 2, . . . of the eigenfunction nodes in z. The dispersion equation for waves belonging to the m-th mode can be written as [21] increases monotonically with k, c and m.
As c reduces (at fixed m and k) and crosses the c = 1 boundary, the wave-flow resonance interaction [13][14][15] comes into play in the critical layer, and the wave loses its stability. Hence, the curve J = J m (k; 1) serves as the upper boundary J = J  Figure 2) in which there appear to be no eigenoscillations belonging to the same mode [21].
The upper and lower boundaries of the m-th instability domain start from the same point of the diagram, J and then increase very slowly ( Figure 2). In what follows, we shall deal with the fundamental mode m = 0 only, and it should be emphasized that, along the lower boundary of the instability domain, its phase velocity differs by only O(ℓ) from U N in a wide range of k, ℓ ≪ k ≪ ℓ −1 , [21]. Therefore, when 0 < J − J * = O(ℓ 2 ), a rather broad spectrum of unstable planar waves is formed with Three-dimensional perturbations can be easily incorporated into this picture employing Squire's theorem [23,8] which states that the problem of evolution (with the time t) of the oblique wave with the wave vector k 3D = (K cos θ, K sin θ, 0) in the flow with velocity and buoyancy profiles U (z) and n(z), and the bulk Richardson number J is equivalent to the problem of evolution, with the time t ′ = t cos θ, of the planar wave with the wave vector k 2D = (K, 0, 0) in the flow with the same U (z) and n(z) but a greater J ′ = J/ cos 2 θ. As a consequence, the complex-valued frequencies of the threeand two-dimensional disturbances are related by Let us now summarize some of the most important properties of Holmboe's instability in flows being studied. First, such a flow is unstable at any stratification level J > 0, and the instability domain contains all the wavelengths (all K > 0). Secondly, so long as J < J * , only oblique (θ ̸ = 0) disturbances are unstable, and they remain the fastest-growing up to J = O(ℓ 3/2 ) (see [21,22]). In the same range, J * < J < ∼ ℓ 3/2 , the flow has a wide three-dimensional spectrum of unstable waves with growth rates γ = O(ℓ) slow depending on k (Figure 3) and nearly equal phase velocities ( Figure 4). As a result, the instability amplifies a wide three-dimensional spectrum of Holmboe waves. And thirdly, critical layers of these waves are at a distance of O(ℓ) from each other and are located within the pycnocline or at its nearest periphery ( Figure 4). When viscosity is small, at the linear stage of development the thickness L of an individual (i.e., belonging to a single wave) CL is equal to the unsteady scale [24], L = L t = O(γ) = O(ℓ), and during the course of further evolution it can only grow. Hence, from the beginning of instability development, individual critical layers of a broad spectrum of waves merge into a common CL, and therefore these waves can efficiently interact during the weakly nonlinear stage of their development. Notice that the pycnocline is, in fact, immersed into such a CL, and they both should be considered as a whole.
In the framework of weakly nonlinear approach, wave interaction is taken into account using the method of successive approximations, and the nonlinear evolution equations (hereafter referred to as NEEs) have the form of expansions in amplitudes of interacting waves (see [25], for example). As long as the amplitudes are small enough, the main contribution to NEEs is generally due to the resonant interaction of triplets of waves whose wave vectors satisfy the triangle condition k 1 = k 2 +k 3 , and frequency detuning, ∆ω = ω 1 − ω 2 − ω 3 , is sufficiently small (this is inherently fulfilled for waves with a common CL). And we shall see that as soon as some amplitude threshold is attained, four-wave (k 1 = k 2 ± k 3 ± k 4 ), five-wave and other higher-order interactions come into play, almost all at the same time rather than one by one. Therefore, it is best to use qualitative methods for describing this stage.
The paper is organized as follows. In Section 2, basic equations are written and their solutions outside CL are found. Section 3 is devoted to constructing the solution inside CL, matching it to the outer solutions, and analysis of nonlinear evolution equations obtained, their solutions and possible scenarios of the perturbation evolution. In Section 4, some numerical solutions of the evolution equations are presented. The results are discussed in Section 5. And some details of derivation of evolution equations are given in Appendix.

Basic equations and outer solution
Let the bulk Richardson number J = ℓ 2J = O(ℓ 2 ), and Reynolds number is so large (Re ≫ ℓ −3 ) that the dissipation is negligible (in numerical simulation of Holmboe instability, Re > ∼ ℓ −4 is usually taken, see [12], for example). We consider the problem of weakly nonlinear development of instability waves belonging to the most unstable, for such a level of stratification, mode m = 0. Because the main nonlinear processes are located within the CL (which is arranged much the same in all flows), we are able to carry out the analysis in general form, rather than for specified profiles of U (z) and ρ 0 (z). The results obtained will be illustrated by numerical calculations for the model flow Since phase velocities of waves involved are close to the mean flow velocity in the pycnocline, U N , in the frame of reference moving downstream with velocity U N , both propagation of waves with phase velocities and their growth in amplitude can be described in terms of slow time τ = Lt where L is the thickness of their common CL. It is assumed that the disturbance evolves from an infinitesimal one due to instability and therefore vanishes as τ → −∞. We shall employ matched asymptotic expansions to construct the solution, i.e., we search for solutions inside CL and in outer regions of the flow, and then, matching them in intermediate domains (L ≪ |z − z N | ≪ 1), we obtain the evolution equations as matching conditions.
Eliminating pressure, one can write the hydrodynamic equations in the Boussinesq approximation [16], in the frame of reference moving with velocity U N , as and ρ are the velocity and density disturbances, and ∆ is Laplacian.
In outer regions of the flow (|z − z N | ≫ L), stratification and nonlinearity can be neglected and one may keep only the linear part of the solution, ] , where summation is over all the waves of the ensemble, and c.c. denotes the complex conjugated term. Since amplitudes of interacting waves can, in general, have (or acquire in the process of evolution) different orders of magnitude, it is convenient to assign them different amplitude parameters 0 < ε n ≪ 1 and to seek the solution for (4) in the form ] .
Outside the CL, the leading-order (O(ε n )) contribution to the n-th wave can be written as w (1) n± = −ik n A n (τ )g ± (z; K n ) where the +/− subscript corresponds to the region over/under z N , K 2 n = k 2 n + q 2 n , and g ± (z; K) satisfy the Rayleigh equation with the boundary conditions Inside the CL, the eigenfunction g(z) of the oscillation belonging to the mode m = 0 changes insignificantly, hence, to the leading order, Coefficients a n+ = a + (K n ) and a n− = a − (K n ) are real, and can be found only by solving the problem (6). It should be emphasized that, for any K n , g + and g − are different solutions of the Rayleigh equation because, according to the Rayleigh theorem, there are no eigenoscillations in homogeneous flows with U ′′ (z) < 0 (exclusive of a continuous spectrum of the so-called Van Kampen -Case modes, see [13,15]). Therefore, the functions g n±(z) ≡ g ± (z; K n ), when analytically continued into the lower half-plane of the complex z, are linearly independent, and their Wronskian Turning back to the vertical velocity, we see that the jump of the w n derivative should be equated with the O(Lε n ) change in the derivative of the W n component of the inner solution. As a result, we arrive at the matching conditions which yield the nonlinear evolution equations.

Inner equations and an outline of their solution
Inside the critical layer, it is convenient to introduce new variables, as well as the time τ ′ = U ′ N τ (in what follows, the prime for τ will be omitted) andR =J/U ′ N 2 . Then, accurate to exponentially small terms (compare with (1)), and (4) take the form (we retain only the terms necessary for further analysis; derivatives in τ, x, y, Z are denoted by corresponding subscripts, f N ≡ f (z N )): ( ] , , ] , Let us begin with a linear part of the inner solution. In line with (5), it is represented as a sum of waves, and each wave as an expansion n , W n , P n , Γ ] .

O(ε n ). Equation (11 a) yields homogeneous equations
n ZZ = 0 and, after matching to the outer solution (7), Other contributions to the inner solution, both linear and nonlinear, obey the equations of the form 134 Weakly Nonlinear Stage of Instability Development in a Sharply Stratified Shear Flow with an Inflection-Free Velocity Profile Substituting (12) into (11 b) and using (13) one obtains and P (1) Note that A Further, we employ (11 c-e) to find the horizontal components of velocity, which, as can be easily verified, are automatically matched to (9). 2. O(L ε n ). In this order, we obtain the equa- The right-hand side is regular in Z = 0. Therefore, its integral in Z does not change if we bypass this point from below. Since the term containing ∂A n /∂τ is regular and decreases as Z −2 (when |Z| → ∞) in the lower half-plane of complex Z its integral in Z vanishes, and ∫ − dZW ) .

3.
O(ε n ℓ 2 /L). In this order, we calculate the correction for stratification, n , and its contribution to (10), Now, we have all the main linear contributions to (10) due to the inner solution. Rearranging them to the lefthand side and keeping nonlinear terms (that have yet to be calculated) on the right, we obtain the equations (see (8)) Before we proceed to calculating W (N ) n , one remark is necessary. We study an ensemble of the most unstable waves. During the linear stage, their amplitudes A n (τ ) obey equations (16) with zero right-hand sides and grow exponentially, A n (τ ) = A n0 exp[−ik n (Z n + i∆ n )τ ], with growth rates γ n = k n ∆ n which satisfy the dispersion equations and are slow-dependent on k n in a broad spectral range (see Figure 3). Hence, to the beginning of a nonlinear stage, a wide spectrum of these waves is formed, with amplitudes of the same order, and it is correct to put the CL thickness L = ℓ and all ε n = ε during not only the linear stage of development but also early nonlinear one. It is easily seen that the nonlinear part of the inner solution can be expanded in ε/L 2 and L, and the leading order of W (N ) n is O(ε 2 /L 2 ) and quadratic in amplitudes, i.e., corresponds to three-wave interactions. Their contribution to the right-hand side of (16) becomes competitive and comes into play as soon as the wave amplitudes achieve the level

Evolution equations and their solutions
Recall that the wave parameters are so defined that the x-components k n of wave vectors as well as frequencies, ω n ≈ U N k n , are positive. Therefore, in the triad k 1 = k 2 + k 3 , the wave k 1 has a maximal frequency. Analysis of the three-wave stage of perturbation development [26] has revealed an interesting feature of the nonlinear wave interaction inside CL. It was found that only 'decay' processes (k 1 − k 3, 2 −→ k 2, 3 ) contribute to (16) whereas 'fusion' processes (k 2 + k 3 −→ k 1 ) do not. As a result, in three-wave interaction the highestfrequency wave plays a catalyzing role [19,20] because it accelerates the growth of two other waves but does not experience their reverse effect and continues to rise with its linear growth rate. And a similar property is inherent to interactions of higher orders. Namely, among m-wave processes k 1 = k 2 ± k 3 · · · ± k m , only those can contribute to the development of the wave k 1 that have at least one minus in their symbolic formulae (a mathematical explanation of this rule is given in Appendix).

Three-wave stage of evolution
As soon as the wave amplitudes achieve the level (18), three-wave interactions come into play, and the result turns out to be dependent on the configuration of wave ensemble [26]. In an isolated triad, the wave k 1 rises with its linear growth rate γ 1 interacting with two Universal Journal of Physics and Application 2(2): 129-141, 2014 135 other waves (the bar denotes a complex conjugate), and parametrically amplifying their growth up to a super-exponential one (∼ exp ). Here If the ensemble contains many waves, the most highfrequency of them does also rise with the linear growth rate and play a catalyzing role. And any of the other waves can be parametrically amplified in some triads and serve as a catalyst in the others, stimulating the waves coupled with it to grow faster. In order to describe them, we extend (19) and construct NEE for some A k (τ ) by inclusion into its right-hand side the sum of contributions of all the triads containing k, where ) .
Solutions of these equations grow much faster, according to an explosive law, Setting k 1 = k + k 2 , we substitute this law into (21) and calculate each of its terms denoting them P 1 , P 2 , and P N , respectively. After some simple algebra we obtain: ) , where Γ(x) is the Euler gamma function, F (a, b; c; x) is the hypergeometric function, and Ψ(a, c; x) is the confluent hypergeometric function [27]. It is easily seen that when τ → τ * from below, |P 2 | ≪ |P 1 |, whereas P 1 and P N have the same dependence on τ if It should be emphasized that, after transition to an explosive growth, the second term on the left-hand side of (21) due to the flow stratification (P 2 ) becomes noncompetitive. This fact is not surprising because the relation of the stratification-induced term in (11 a) to the following one is of the order O(ℓ 2 /L 2 ) and L increases with τ , Now let us impose the (commonly accepted in numerical simulation) condition of disturbance periodicity, with periods 2π/k in x and 2π/q in y. Then the ensemble should consist of waves with wave vectors of the form (mk, pq, 0), with integer 0 < m ≤ M and p. In particular (hereafter, zero z-components will be omitted), the wave vector of the highest-frequency wave k M = (M k, p M q) (M ≥ 3 because M = 2 corresponds to isolated harmonic-subharmonic triad(s) that grows super-exponentially rather than explosively). In what follows, such an ensemble will be called the M -level ensemble, and the levels will be numbered by m. Waves belonging to the same level have nearly the same frequencies and differ in p (i.e., in the y-component of the wave vector). We assume that, asymptotically, the wave amplitudes at all levels grow explosively with indexes α m depending on m only, and formally assign α M = 0 to the catalyzing wave. It should be also mentioned that the first condition (23) is sufficient but not necessary. Indeed, the growth rate of P N is obviously determined by the maximal sum of indexes of interacting waves, α N = [max k2 (α k1 + α k2 ) − 3], therefore, P 1 and P N grow at the same rate when α k = α N . Applying this condition to our ensemble yields Thus, the fastest growing waves are the first-level waves, and growth index decreases with m (with the step ∆α) to zero at m = M . Numerical calculations of the threewave stage of development in ensembles with M = 3, 4 and 5 [26] have demonstrated that wave amplitudes grow explosively, with growth indexes (25). In particular, in ensembles with M = 3 the growth index of the firstlevel waves is twice as large as that on the second level, α 1 = 6 and α 2 = 3.

Post-three-wave stage: a qualitative analysis of interactions
Thus, in sufficiently wide ensembles of waves, threewave interaction leads to an explosive growth of amplitudes. The CL thickness (24) increases in the process (still remaining small compared to unity), and the scaling (18) should be replaced by relations Using them, one can estimate the contributions of higher-order processes into the right-hand side of (16) for the wave belonging to the m-th level. But first recall that the contribution to W (N ) due to the interaction of waves with amplitudes O(ε 1 ) and O(ε 2 ) has an amplitude O(ε 1 ε 2 /L 2 ) (see (11 a)), and that only those processes contribute to NEE which have at least one minus in their symbolic formulae. In three-wave interaction (the levels of interacting waves are numbered by m i ), the maximum contribution to (16) is as expected. Further, it is easily seen that only two types of fourwave interactions can contribute to the evolution equations for waves of the m-th level, m = m 1 + m 2 − m 3 and m = m 1 − m 2 − m 3 . In the process of the first type, to the same m 2 and m 3 a wave with a smaller m 1 corresponds, and therefore the contribution of such a process is greater. A simple estimate, shows that this contribution becomes of the order of unity when ) , or, expressing ε m in terms of L, Assuming that these amplitudes are reached, we are able to evaluate the contributions of further orders. As earlier, we can conclude that the maximal n-wave contribution is due to the interactions of the form k = k 1 + k 2 + · · · + k n−2 − k n−1 . Taking the x-component of this equality, we find and, using (28), obtain max εm 1 εm 2 ...εm n−1 εmL 2n−3 = max L 2(n−1)+(m 1 +m 2 +···+m n−1 )/2 Thus, upon reaching the amplitudes (27), many interactions of various orders come simultaneously into play. Namely, all the interaction from the three-to the (m + 3)-wave inclusive begin to contribute to the development of each wave belonging to the m-th level. In general, the exception are only the most high-frequency waves, i.e., waves of the M -th level, because their evolution is not affected by three-wave interaction. According to (30), the interactions of different orders no longer compete with each other and equally affect the development of disturbances. Therefore, the scaling (28) is maintained until the end of the weakly nonlinear stage of development, i.e., up to L = O(1) when all the wave amplitudes on all levels also become of the order of unity.
It is interesting to note that (27) is a threshold level for isolated triads of interacting waves as well (for them, M = 2). As soon as the amplitudes of two lowerfrequency (or the 'first-level') waves, k 2 and k 3 , which increase super-exponentially become O(ℓ 5/2 ), they begin to affect the development of the third (higher-frequency, m = M = 2) wave k 1 by four-(k 1 = k 1 + k 2,3 − k 2,3 ) and five-wave (k 1 = k 2 + k 3 + k 2,3 − k 2,3 ) interactions, and they all begin to grow explosively, with indexes 3 and 5/2, see (28).

Post-three-wave stage: the structure of evolution equations
Since a large number of nonlinear interactions of diverse orders take part in the perturbation development, the deriving the evolution equations is extremely timeconsuming, and the equations themselves are highly complicated. Based on qualitative considerations, we try to recognize what the nonlinear interactions should be taken into account and what the structure of NEEs should be. In the process, we shall appeal to Appendix where a rather detailed derivation of the three-wave contribution is given with appropriate comments.
After transiting to an explosive growth, a distinctive hierarchy of amplitudes is set up (see (26) and (28)). Namely, the higher is the frequency of the wave (i.e., the level number m), the less is its amplitude, and only to the end of a weakly nonlinear stage all the amplitudes reach the same order of magnitude. Therefore, only that n-wave interaction is efficient which has one minus and m n−1 = 1 on the right-hand side of its formula (29), and the remaining processes are uncompetitive. Moreover, at an explosive stage the linear contributions to NEEs due to the stratification become uncompetitive as well (see the text before (24)).
Recall, however, that there are also the preceding stages of development, namely, the linear stage and the transition from it to an explosive growth. For the entire evolution to be described correctly, it is necessary to retain in (16) all the linear and three-wave contributions, and only those contributions of higher-order interactions that are of the order of unity.
Let us turn to the structure of such contributions. We start with a scheme of calculation of the inner solution. The iteration process begins with Z-independent leading-order solutions (12). Each further iteration obeys the equation of the form L n f (τ, Z) = R(τ, Z) and is given by an integral over the time delay s (13) which depends on all the history of perturbation development. Since R(τ, Z) is expressed by an integral of the same sort, f (τ, Z) has the form of a multiple integral (or their sum) from 0 to +∞ over all the variables, with Z multiplied by a linear combination of s, s 1 , . . . , s p in the argument of exponential function. It is easily seen that differentiation with respect to Z is reduced to the multiplication of F by (−i) and the same linear combination. Hence, as one can see from (11 a), in each following order of interaction, W (N ) ZZ acquires one integration and its integrand is multiplied by a linear combination of s i and an additional amplitude function. Note that since each amplitude has its own time delay, the number of integrations cannot be less than the number of amplitude functions. As a result, a n-wave contribution to W (N ) ZZ is a homogeneous function of degree (n − 1) with respect to amplitudes and of degree 2(n − 1) with respect to s i (see (A1) and (A2) for comparison).
Let us take a wave of the m-th level and estimate the contribution W ′′ mn to W (N ) ZZ made by an efficient n-wave interaction (29) with m n−1 = 1. It is a sum of terms that differ in structure but all are homogeneous functions of the same degrees of amplitudes and s i . A single term is a (r + n − 1)-fold integral (0 ≤ r < n − 1) of the form , s b , . . . , s r ; s, s 1 , . . . , s n−2 ) where µ = ms , the kernel K is a homogeneous function of degree (n − 1 − r) of its arguments, each S j is a sum of some of the delays s i (i = a, b, . . . , r, 1, 2, . . . , n − 2), possibly empty, and all S j are different (see (A1) and (A2)). Such a structure is typical for all higher-order interactions, from four-to (m − 3)-wave.
Substituting W ′′ mn into (16), integrating it over Z, and using the relation Let now find the τ -dependence of this contribution to NEE when the amplitudes grow explosively, A i (τ ) ∼ (τ * − τ ) −αm i with α m = 2 + m/2, see (28). Introducing Hence, all the higher-order nonlinearities exert a joint action on the disturbance rather than compete with each other.

Numerical study of evolution
In order to illustrate the above analysis, we have calculated a weakly nonlinear evolution of 'minimal' threelevel (M = 3) ensemble k 0 = (3k, q), k 11 = (k, −3q), k 12 = (k, 2q), containing two by two waves at the first and the second levels, and one wave at the third. Drawing up the scheme of interactions for this ensemble, we include in it all the three-wave interactions with a non-zero contributions, and only efficient higher-order processes. Also, it is taken into account that the interactions involving only the waves with collinear wave vectors are uncompetitive (for discussion, see [28,29], for example), and that k 12 and k 22  It can be seen that even in such a minimal ensemble the number of interactions that should be taken into account is immense (compare with the three-wave stage when each wave takes part only in one or two interactions, see the scheme and [26]).
We calculate all the coefficients of interaction (they are not written here because of their large number and awkwardness) and find a numerical solution of the evolution equations for k = 0.3 and q = 0.5. The numerical procedure is the same as described in [26]. Figure 5 shows the position of the ensemble waves inside the instability domain, and Figure 6 demonstrates the τ -dependence of the wave amplitudes.
In view of the above analysis, as soon as amplitudes attain the O(ℓ 3 ) level, the linear stage of development (corresponding to rectilinear parts of amplitude traces) changes to the three-wave one, and both the first-and the second-level waves begin to grow explosively, as (τ A − τ ) −αm with α 1 =6 and α 2 =3, whereas the thirdlevel wave continues to grow exponentially and interacts with other waves as a catalyst only. And when the wave amplitudes at the levels with m = 1 and 2 become as high as (27), all the waves begin to increase explosively as (τ B − τ ) −(2+m/2 , i.e., with α 1 = 5/2, α 2 = 3 and α 3 = 7/2. On the Figure 7, the explosive stages are shown in more detail, and straight lines with inclinations corresponding to α mentioned above are drawn for comparison. One can see that the numerical results are in a good agreement with the predictions. In conclusion, it should be noted that the three-wave and post-threewave stages have different times of 'explosion', τ A and τ B respectively. In addition, we have calculated the evolution of isolated wave triads in two versions, with a frequency ratio of 1:2:3 (waves with k 12 , k 21 and k 0 , see (32)) and with an 'oblique' subharmonic (k 1 = (0.9, 0.5), k 2 = (0.45, −0.3), k 3 = (0.45, 0.8)). The results obtained are in a good agreement with the theory as well.

Discussion
The flows of the class under study (i.e., sharply stratified and having no inflection points on the velocity profile) have an outstanding distinction (as compared with a homogeneous boundary layer, for example) giving a better insight into the mechanisms operating on the way from linear instability to a strongly nonlinear stage of perturbation development, i.e., into the mechanisms of transition to turbulence. The distinction is that the instability waves in those flows belonging to the funda- mental mode (m = 0) have very close phase velocities (and therefore, a common critical layer) in a large part of their spectrum. Moreover, the pycnocline turns out to be immersed into the unsteady critical layer with a rapidly changing flow structure. As a result, both the dispersion and nonlinear terms of evolution equations depend on the whole history of the disturbance development, hence, evolution equations have the form of a system of integral equations, see (16) and (19). We have analyzed the overall stage of a weakly nonlinear evolution and have found that in the course of its second, post-three-wave part, a large number of wave interactions should be taken into consideration. As a result, nonlinear evolution equations become too cumbersome, and their deriving seems to be an almost impossible task. Therefore, in this study we relied on a qualitative analysis of the problem and, owing to it, were able to find the structure of the evolution equations and to analyze the behavior of their solutions. And only for a few very small wave ensembles we have derived all the nonlinear interaction terms and numerically calculated the evolution of wave amplitudes (see Section 4).
By analogy with the acoustics, it could be expected that in the ensemble of waves with very close phase velocities (i.e., weak dispersion), the rapid growth of high-frequency harmonics accompanied by formation of fronts and structures similar to weak shock waves [30,31] should take place. However, during the initial (threewave) stage of nonlinear development the scenario is exactly the opposite. Namely, the fastest-growing is a lowfrequency component of the spectrum (waves with small m, see (25) and (26)), and the disturbance becomes even more smooth. The reason for this is in a catalytic character of three-wave interactions within the critical layer where the wave fusion (k 2 +k 3 −→ k 1 ) is forbidden. And only on the next (post-three-wave) stage, when a lot of interactions of diverse orders come into play, the growth of high-frequency part of the spectrum is accelerated. Nevertheless, only to the end of weakly nonlinear stage all the amplitudes become of the same order. Note that over the weakly nonlinear stage all the wave amplitudes grow explosively (with changing their growth indexes in the course of transition from three-wave stage to postthree-wave one).
In conclusion, let us dwell on one more issue. One of the most significant achievements of weakly nonlinear theory is calculation of the spectra of weak turbulence, in particular, the spectra of wind waves [32], which are in a reasonable agreement with observations. It would be interesting to calculate such spectra for Holmboe waves and compare them with Kolmogorov spectrum for developed turbulence. Unfortunately, as the above analysis demonstrates, Holmboe waves grow too rapidly (explosively) and have no sufficient time for phase mixing. They 'remember' the initial conditions up to the end of the weakly nonlinear stage of their development, and turbulent spectrum can not form during this stage. Derivation of the three-wave NEE For three waves with k 1 = k 2 + k 3 , using (12), (15) and notations (20) we derive from (11 a) ] , 3 ) + 2 ) + ] .
It is easily seen that the equations for W can be obtained from each other by exchange of the indexes 2 and 3. Therefore we write out only the result of integration for the first and the second equations keeping in mind that k 1 = k 2 + k 3 : (A2) Note that the right-hand sides of (A1) and (A2) are homogeneous functions of degree 2 with respect to amplitudes and of degree 4 with respect to s i . In arguments of exponential functions, the expressions in square brackets are sums of ±k i S i where S i is the delay of A i , and only items due to complex conjugated amplitudes are negative (compare with (29) and (31)). Those exponential functions are underlined in which the expression in square brackets vanishes if and only if the s i appearing in it are all zero. On integrating over Z, the terms proportional to such exponents vanish. In particular, the integral of W which is homogeneous function of degree 2 with respect to amplitudes and only of degree 3 with respect to s i . Finally, let us call attention to the fact that the term containing complex conjugated amplitudes need not make a nonzero contribution to (16). The second item on the right-hand side of (A2) is a good illustration.