Evaluation of a spectral line width for the Phillips spectrum by means of numerical simulation

The work aims to check one of the assumptions under which the kinetic equation for water waves was derived in order to understand whether it can be applied to the situations described by the Phillips spectrum. We evaluate a spectral line width of the spectrum from the simulations in the framework of primordial dynamical equations at different levels of nonlinearity in the system, corresponding to the weakly turbulent Kolmogorov–Zakharov spectra ω, Phillips spectra ω, and intermediate cases. The original motivation of the work was to check one of the assumptions under which the kinetic equation for water waves was derived in order to understand whether it can be applied to the Phillips spectrum. It is shown that, even in the case of relatively high average steepness, when the Phillips spectrum is present in the system, the spectral lines are still very narrow, at least in the region of the direct cascade spectrum. It allows us to state that, even in the case of the Phillips spectrum, one of the assumptions used for the derivation of the Hasselmann kinetic equation is still valid, at least in the case of moderate whitecapping.


Introduction
Two really seminal papers were published in the area of physical oceanography about half a century ago: the article by Phillips (1958) and the work of Hasselmann (1962).Both of them concentrated on the same problem: what is going on with ocean waves growing under the influence of wind?O. Phillips suggested that this growth is arrested by wave breaking, in other words, by the formation of "white caps" (or "white horses").Wave breaking is the main mechanism of energy dissipation.In "white caps", mechanical energy of waves transforms to turbulence on small scales, then to heat.This is a strongly nonlinear phenomenon that cannot be studied by perturbative methods.An analytical theory of wave breaking has not been developed yet.The Phillips' assumption on dominance of wave-breaking effects allowed him to offer, on the basis of dimensional considerations, a universal ω −5 energy spectrum of wind waves (Phillips' spectrum), whereas a more elaborate explanation of the underlying physics was put forward by Kuznetsov (2004).
K. Hasselmann developed a completely different theory.He noticed that a typical ensemble of ocean waves possesses a small parameter: the average steepness (slope) of the surface µ.One can define it, for instance, as where η is the shape of the surface.This definition has the clearest geometrical meaning.Another definition estimates it as the average amplitude A of the waves multiplied by the wave vector of the spectral peak k p .This definition is easy to use for experimental observations.Later, we shall consider another definition, which involves the spectral distribution function.All these definitions give close values of the average steepness.Typically, µ 0.1.This fact allows one to assume that the wind-driven sea is an ensemble of weakly interacting waves, which can be described statistically by the Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
A. O. Korotkevich and V. E. Zakharov: Evaluation of a spectral line width for the Phillips spectrum use of expansion in powers of average steepness µ.This is a tedious procedure because one has to expand up to order µ 4 .However, K. Hasselmann coped with this hard work and derived his famous Hasselmann kinetic equation for a wave energy spectrum.Then, it was found that the Hasselmann equation has exact stationary solutions (Kolmogorov-Zakharov or KZ spectra) (Zakharov and Filonenko, 1967a, b;Zakharov and Zaslavskii, 1982;Zakharov et al., 1992) decaying at a highfrequency range as ω −4 , i.e., slower than the Phillips spectrum.A question arises as to which theory is closer to reality.The answer is given by experiment.
Numerous measurements (e.g., Hwang et al., 2000;Donelan et al., 1985) performed in large lakes and the ocean demonstrated that the spectrum is the weak-turbulent KZ one E(ω) ∼ ω −4 , at least for the range of scales ω p < ω < 4ω p ; here, ω p is the frequency of the spectral peak.At the same time, the Phillips spectrum was observed for rough seas for the tail of the spectrum.This fact has the following explanation.Wave turbulence in the wind-driven sea is a mixture of weak and strong turbulence.This is a question of phase correlation.Weakly nonlinear interaction provides correlations of waves phases on large scales -hundreds of characteristic wavelengths.At the same time, the wave breaking and whitecapping is the phenomenon localized in space.As a result, looking at the wind-driven sea, one observes the formation of short-lived localized wave breaking events embedded into a homogeneous weakly nonlinear background (Nazarenko et al., 2010).
A similar situation, the coexistence of weak turbulence and localized coherent structures, is typical of wave turbulence in general.It takes place, for instance, in nonlinear optics (Dyachenko et al., 1992), where optical turbulence coexists with self-focusing or in isotropic plasma where weak Langmuir turbulence coexists with Langmuir collapses (Zakharov et al., 1992).Similar phenomena were observed recently in numerical simulations for gravity waves (Zakharov et al., 2007).In all these situations, the usual course of action is to augment corresponding kinetic equations of weak turbulence by the introduction of an additional term, describing dissipation of energy in the coherent structures (corresponding to whitecapping onsets and wave breaking) (Zakharov et al., 2009).The question arises whether it is still reasonable to use a kinetic equation as an adequate model of wave interaction in the presence of coherent structures.Maybe, the first attempt to answer this question was the work by Hasselmann (1974) on spectral dissipation due to whitecapping.In this work, the concept of a "weakly nonlinear on average" sea was introduced, which justifies the weakly nonlinear description, resulting in the Hasselmann equation.This hypothesis was not confirmed numerically because the required numerical simulation of primordial dynamical equations was not feasible at that time.
In this paper, we show that at least some of the assumptions under which the Hasselmann kinetic equation was de-rived still hold in the presence of a number of whitecapping events.To do this, we perform massive numerical simulations of primordial weakly nonlinear dynamical equations (Hamiltonian equations with a Hamiltonian truncated at fourwave interaction terms) describing potential flow of Euler equations of an ideal fluid with a free surface.It should be noted that we only mimic energy dissipation due to whitecapping through the regularization of the equations by relatively close dissipation regions (approximately one decade from the pumping scale).As a result, although the whitecapping cannot be described rigorously in the weakly nonlinear model, one can dissipate part of the energy of the wave involved in the whitecapping by the proposed mechanism.The statement that we dissipate the right part of the energy is a hypothesis, although we give some explanations showing that it looks like quite a reasonable claim.
In our numerical experiments, we study Fourier spectra of space-time correlation functions of normal canonical variables describing the surface dynamics.Our goal was to determine the "shape of the line", i.e., the frequency spectra of certain spatial Fourier harmonics.According to the theory of weak turbulence, these spectra must be narrow and concentrated near the frequency given by the dispersion relation for waves of small amplitude.We obtained the following important result: in a broad range of wave numbers, the spectra remain narrow even in the presence of whitecapping dissipation events, when the spatial spectrum obeys the Phillips law.This justifies the construction of an appropriately adjusted dissipation term that can be used thereafter in wave forecasting.

Basic model
We solve numerically weakly nonlinear Euler equations for dynamics of incompressible deep fluids with free surfaces in the presence of gravity by the pseudo-spectral code described in Korotkevich et al. (2012).The code was verified in our previous papers (Dyachenko et al., 2003a(Dyachenko et al., , b, 2004;;Zakharov et al., 2005Zakharov et al., , 2007;;Korotkevich et al., 2008;Korotkevich, 2008Korotkevich, , 2012Korotkevich, , 2013)).The equations are written for surface elevation η(x, y, t) and hydrodynamic velocity potential on the surface ψ(x, y, t).Equations are the result of weakly nonlinear expansion of the Hamiltonian up to the fourth-order terms in steepness (Zakharov et al., 1992) Here, the dot means the time derivative, = ∇ 2 the Laplace operator in the r = (x, y) T plane, k = √ − , F −1 is an in-verse Fourier transform, γ k is a dissipation rate depending on the scale (according to recent work by Dias et al. (2008), it has to be included in both equations), which corresponds to viscosity on small scales and, if needed, "artificial" damping on large scales.P k is the driving term that simulates pumping on relatively large scales.Here and further on, we use the following formula for coefficients of Fourier series The specific form of pumping is not important in our simulations.We have chosen it in the following form: (3) Here, k p1 = 28, k p2 = 32 and F 0 = 1.5 × 10 −5 ; R k (t) was a uniformly distributed random number in the interval (0, 2π] for each k and t.The initial condition was a low-amplitude noise in all harmonics.The time step was t = 6.7 × 10 −4 .The dissipation function is k was the same in all experiments: where k d = 256 and γ 0 = 0.97 × 10 2 .γ (2) k is the dissipation concentrated in small wave numbers.It was zero only in one of three experiments.It is described separately for every experiment below.
The linear dispersion relation for gravity waves was ω k = √ g k, with g the gravity acceleration, which was taken as g = 1.Equation ( 2) was solved in a 2π × 2π periodic box.The number of modes was 1024 × 1024.

Derivation of contribution to steepness from different scales
Let us consider explicitly discrete Fourier transforms in the case of the 2π × 2π periodic box: Here, we took into account the fact that, for the 2π × 2π region, k = 1.
The same for ψ r and ψ k .Normal canonical variables are Inside, the box spectra were practically isotropic; thus, one can put For average squared steepness, we obtain Using Eq. ( 6), one can get If we consider cases of random phases (as in our simulation), in the isotropic case, where • φ means averaging over the azimuthal angle Using Eq. ( 7), we can simplify the expression for µ 2 : For isotropic cases, it is natural to introduce a spectrum averaged by an angle (this is what we usually measure in our numerical experiments): Then, we have www.nonlin-processes-geophys.net/22/325/2015/ Nonlin.Processes Geophys., 22, 325-335, 2015 For surface gravity waves over the fluid of infinite depth, ω k = √ gk.In our simulations, g = 1, so from Eq. ( 9), we get We would like to emphasize that, for the isotropic case definition, Eq. ( 10) coincides with Eq. ( 1).
In order to understand which scales are important for the steepness, we can introduce µ 2 (k) as follows: Essentially, we introduced a steepness µ(k) accumulated by the scales with wave numbers less than k as a square root of Eq. ( 11).
Here, we use angle-averaged n k (Eq.8), and summation is replaced by integration over k = |k|.Average steepness can be obtained as a limit: This definition gives a value of average steepness that differs from the geometrical definition µ = |∇η| 2 ; in our experiments, the discrepancy does not exceed 0.005, which is in most cases just a few percent of the characteristic average steepness.This discrepancy arises from the fact that, close to the origin of the k plane, we have a small number of harmonics for the angle averaging.As a result, terms in Eq. ( 7) do not vanish completely, because the summation over a coarse discrete grid of wave numbers is too crude an approximation of the integration.
It is worth noting that the usual definition of the steepness through the average squared amplitude and the wave number of the spectral peak takes into account mostly regions of the spectral peak, while our computations of function µ(k) show that the contribution from the tail of the spectrum is a significant part of total steepness µ and cannot be neglected.
For the spatial harmonic with wave vector k, we defined the time Fourier transform as follows: Here, T = 2π/ω k is the period of the chosen wave.Thereafter, I (k, ω) = |a(k, ω)| 2 , which is a spatio-temporal spectrum of the wave field.

Description of experiments
We performed three series of experiments, choosing different functions γ (2) k describing damping in the area of small wave numbers.
The picture we observed is typical of wave turbulence and was described theoretically in detail by Zakharov et al. (1992).Firstly, we get fast growth of waves corresponding to the pumping scale.Then, a direct cascade of energy results in the relatively fast formation of the spectrum tail from the pumping region to the high wave numbers until it reaches the dissipation scale k d .The wave turbulence theory predicts the KZ tail for the direct cascade . We observed both spectra.
There is a simultaneous, although much slower, process of formation of the inverse cascade, corresponding to the flux of wave action from the pumping scale into the larger scales (smaller wave numbers).After very slow development of the inverse cascade spectrum, formation of a strong long wave background starts.This strong long wave background is often called the condensate (the name comes from some similarity to the phenomenon of Bose-Einstein condensation).Formation of the condensate is due to the fact that, at large scales (small wave numbers), the discrete grid of wave vectors becomes too coarse for resonant interaction of waves.In this range of scales, we switch from normal weak turbulence through the mesoscopic regime as described in Zakharov et al. (2005) into regimes without nonlinear resonant interaction.As a result, we have an influx of waves' action from small scales to the large scales due to the inverse cascade, which is arrested by the discreteness of the wave number grid at large scales.We get an accumulation of waves on the border of the scales where at least mesoscopic wave turbulence is still possible.These accumulated waves form the intense long wave background, which we call the condensate.It should be emphasized that a discrete grid of wave vectors appears not only in numerical simulations, but also in all wave tank experiments or other finite-size systems.The formation of the condensate in large lakes or the sea can be neglected due to the fact that downshift of the spectral maximum becomes slower the further it is shifted in large scales.As a result, the free path time of the wave is not enough to reach the scale where resolution of the grid of wave numbers becomes significant.

Without condensate and inverse cascade
In the first series of experiments, we assumed that γ (2) k was linearly growing in small wave numbers k < 28: In this case, the inverse cascade of wave action was completely suppressed.We observed the formation of the direct cascade, reasonably described by the standard KZ spectrum |a k | 2 β 1 k −4 , which corresponds to energy spectrum E(ω) ∼ ω −4 , with β 1 1.2 × 10 −4 .This spectrum was observed in the range of scales 32 < k < 150.The spectra in the linear and logarithmic scales are presented in Fig. 1.Here and further on we consider angle-averaged spectra of direct and inverse cascades.This averaging can be considered to be a cheaper replacement for ensemble averaging in the case of isotropic spectrum with respect to angle situation.The steepness grows in this range of scales as shown in Fig. 2. At k max 200, steepness reaches its saturation level µ max 0.104.The shapes of spectral lines at k equal to 50, 100, 150, 200, and 250 are represented in Figs.3-7.One can see that the spectral lines are narrow.Only in areas of large k 250 do we observe bounded (slave) harmonics that do not obey the linear dispersion relation.

Without the condensate, an inverse cascade is present
In the second series of experiments, the low-frequency damping was presented only in small wave numbers k < 10, where the dissipation rate was a constant γ periment, we observed the formation of both direct and inverse cascades.The dynamic range for inverse cascades was too short for accurate evaluation of the power of the spectrum; the only statement that can be made is that it is relatively close to the predicted |a k | 2 ∼ k −23/6 spectrum, with the power exponent α in k −α ranging from 3.1 to 4.0.In the area of direct cascades, we observed a slightly steeper spectrum than the KZ-spectrum slope.In the area 32 < k < 150, the spectrum can be approximated by |a k | 2 ∼ k −4.2 .This deviation from the pure KZ slope k −4 can be explained by the influence of wave breaking and whitecapping effects (Korotkevich, 2012).The spectra are presented in Fig. 8.The steepness in this case reaches the level µ max 0.130 (see Fig. 9).The spectral lines in this series of experiments were still narrow for k = 50, 100, 150, and 200.In the area of significant damping k = 250, we observed intensive formation of slave (bond) harmonics (see .We would like to mention that, although for harmonics around k = 250, dissipation is not formally present (it starts from k d = 256), in reality, dissipation manifests itself in the inertial interval through the process of multiple harmonics generation.If we consider a second harmonic generated from a harmonic around k d = 250, it will be in the region of very strong dissipation and, as a result, energy will be drained from the harmonic formally free from dissipation.In other words, additional dissipation appears in the inertial interval.It has to be noticed that, for sharper waves, this nonlinear process is more efficient, which means the drain of energy through multiple harmonics generation is stronger.This explains the noisy character of the spectral lines' plot, which is even more prominent in the next experiment due to stronger nonlinearity.
It should be noted that, when we say "sharper waves" above, it does not mean that we violate the weakly nonlinear approximation under which our model equations were derived.When the wave starts to sharpen, the spectrum widens and reaches the dissipation scale.Because, in our model dissipation, the scale is relatively close to the pumping (only approximately one decade away), this dissipation regularizes the wave by draining the energy from it through the mechanism described above and, as a result, making the wave less sharp.So, even when waves start to sharpen after some threshold, the spectrum widens enough and the process stops, preventing violation of the assumptions used for the derivation of primordial Eq. (2).

Both the condensate and inverse cascade are present
In the last series of experiments, we completely eliminated low-frequency dissipation γ (2) k = 0.In this case, we observed the formation of an intensive inverse cascade leading to the creation of the condensate at k 5 − 4.
The spectrum in the area of the inverse cascade is the same as in the previous series of experiments, while in the area of the direct cascade, we observed the Phillips spectrum Korotkevich, 2008).The spectra are presented in Fig. 15.In order to check that the spectrum slope is not changing any more, long time calculation was performed on a smaller grid, which showed motion of the condensate position by one wave-number grid step without any notable influence on the inverse cascade region (dashed line in Fig. 15).The average steepness was significantly higher than in the previous experiments, reaching its maximum value of µ max 0.142 (see Fig. 16).We would like to stress that the difference from the previous case, which looks quite small in absolute values (0.130 and 0.142), is really quite strong, because in this range the dissipation rate depends strongly on the average steepness (Zakharov et al., 2009); the probability of whitecapping grows pretty quickly with average steepness (Banner et al., 2000).It is important to notice that the steepness of the condensate was quite moderate (µ 0.06; it is worth noting that average steepness gets a significant contribution from small scales, quite far from the spectral peak, which means that the often used definition of the average steepness through the product of mean amplitude and wave number of the spectral peak can deviate significantly from the geometrical definition through the average slope of the surface).However, modulation of the short waves by long waves caused whitecapping that resulted in the Phillips spectrum (the process described qualitatively in Korotkevich (2008Korotkevich ( , 2012))).The physical picture is quite simple: shorter waves propagating over a strong long wave background (condensate) that moves with acceleration tend to increase the steepness of the front slope, which can result in the onset of whitecapping, causing additional dissipation.One could recall that our governing primordial dynamical Eq. ( 2) was derived in the weakly nonlinear approximation, while foam formation in the case of whitecapping is a strongly nonlinear phenomenon.It could be mentioned that the description of wave breaking in all details would give us a picture of how the energy was transferred from large to small scales and dissipated.At the same time, if we need only to dissipate the wave energy transferred to small scales during foam formation, one can use the dissipation already available in the scheme for robustness.Although we shall not be able to describe the process of dissipation of energy in all details, most of the energy, which was transferred into high wave numbers when one of the slopes of the wave became too steep, will be dissipated.Here is a simple analogy: let us suppose we have pile of stones on a cliff and throw them down one by one; in order to understand how many stones fall down, all we need to do is to calculate how many of them have been thrown; we do not need to describe the rotation of every stone in all details while it is falling down.As a result, we mimic whitecapping dissipation, staying in the framework of a weakly nonlinear Hamiltonian expansion with all the advantages of an accurate description of the nonlinear wave interaction, without fully nonlinear simulation of the wave breaking, which is challenging even in the case of 2-D hydrodynamics (Chalikov and Babanin, 2012).Similar stimulation of whitecapping dissipation can be caused by interaction of short waves with inhomogeneous background currents or strong long-scale swell in the ocean.
The most interesting and important question we focus upon is about the shape of spectral lines in the area of the Phillips spectrum.Samples of spectral lines are presented in Figs.17-21.For this experiment, we analyzed the longest time series that resulted in higher-frequency resolution.After obtaining the frequency spectrum, we used moving averaging in order to get rid of noise.
One can see that, in the most interesting area 32 < k < 150, the spectral lines are still narrow, while essentially "contaminated" by the slave harmonics.The maxima of the spectral peaks are shifted to the high-frequency area according to the theoretically predicted nonlinear frequency shift.In the area of shorter waves k > 200, the spectrum is a chaotic mixture of leading and slave harmonics.A detailed picture characterizing how the dispersion relation changes in the presence of the condensate can be found in Korotkevich (2013).

Conclusions
We believe that our experiments make possible a "Marriage of Heaven and Hell" in the spirit of William Blake.Both outstanding scientists -O.Phillips and K. Hasselmann -were right.If the local steepness µ is small (µ ≤ 0.1), the Hasselmann kinetic equation is valid without any augmentation by additional dissipation terms, at least for not too broad a spectrum, similar to what we have in the experiments.For µ noticeably higher than 0.1, in the area of very short waves, the kinetic equation is not applicable in its "pure" form.It is impossible to separate "free" and "slave" harmonics in this area.This part of the ocean spectrum cannot be described analytically by the use of perturbative methods.Nevertheless, the Phillips spectrum in this area is still applicable.The only theoretical reason for this statement is dimensional considerations supported by experimental data.At the same time, there is also a "grey area", which is especially interesting be- cause it generates a significant part of the steepness where whitecapping has equal footing with weakly nonlinear resonant wave interaction.As a result of the balance of nonlinear flux and dissipation, the spectrum slope can change (Newell and Zakharov, 2008), as has been demonstrated in our simulations.
We expect that spectra in this range could be described by an "augmented" Hasselmann equation, including an additional term describing dissipation energy due to wave breaking.Similar additional "dissipation terms" S diss are widely used in well-developed operational wave forecasting models.But, they are introduced "out of the blue" and are not supported either by theoretical considerations or, in many cases, by experimental observations.Moreover, it is not clear a priori that one can use the Hasselmann kinetic equation in a situation where wave-breaking events are frequent enough.We hope that our experiments showed that this is likely possible.A necessary condition of the applicability of the Hasselmann equation is narrowness of the spectral line (we should stress that, in this paper, the "dynamical" effects, like "gusty" pumping, described in a recent work by Annenkov and Shrira (2009) and confirmed numerically in Annenkov and Shrira (2011), are not considered).In the present article, we assert that, in the "grey area" with a frequency range in one half of the decade, the frequency spectra of harmonics are still narrow lines.It means that the Hasselmann equation is applicable there in the sense that it is a good approximation for the description of nonlinear wave interaction.Of course, it must be augmented by a proper dissipative term.The existing and widely used dissipative terms can be improved.One possible and hopefully plausible variant was recently suggested by Zakharov et al. (2012), another by Zakharov and Badulin (2015).The dependence of the dissipation term on the average steepness was recently measured directly from the numerical experiment; the preliminary results can be found in Zakharov et al. (2009).In our next paper, we shall analyze what information about possible shapes of the dissipative term can be extracted from massive numerical experiments.

Figure 1 .Figure 2 .
Figure 1.Angle-averaged spectrum with both the condensate and inverse cascade suppressed.Slopes for the KZ spectrum and Phillips spectrum are shown for reference by dashed and dotted lines respectively.Main figure: logarithmic scale; inset: linear scale.

Figure 9 .Figure 10 .
Figure 9. Steepness as a function of absolute value of wave number.Only the condensate was suppressed.

Figure 15 .Figure 16 .
Figure 15.Angle-averaged spectrum with the condensate.Two grid sizes.Slopes for the KZ spectrum and Phillips spectrum are shown for a reference by dashed and dotted lines respectively.Main figure: logarithmic scale; inset: linear scale.