Ocean swell within the kinetic equation for water waves

Effects of wave-wave interactions on ocean swell are studied. Results of extensive simulations of swell evolution within the duration-limited setup for the kinetic Hasselmann equation at long times up to $10^6$ seconds are presented. Basic solutions of the theory of weak turbulence, the so-called Kolmogorov-Zakharov solutions, are shown to be relevant to the results of the simulations. Features of self-similarity of wave spectra are detailed and their impact on methods of ocean swell monitoring are discussed. Essential drop of wave energy (wave height) due to wave-wave interactions is found to be pronounced at initial stages of swell evolution (of order of 1000 km for typical parameters of the ocean swell). At longer times wave-wave interactions are responsible for a universal angular distribution of wave spectra in a wide range of initial conditions.


Physical models of ocean swell
Ocean swell is an important constituent of the field of surface gravity waves in the sea and, more generally, of the sea environment as a whole. Swell is usually defined as a fraction of wave field that does not depend (or depends slightly) on local wind. Being generated in relatively small stormy areas these waves can propagate at long distances of many thousand miles, thus, influencing vast ocean stretches. For example, swell from Roaring Forties in the Southern Ocean can traverse the Pacifica and reach distant shores of California and Kamchatka. Predicting of swell as a part of surface wave forecast remains a burning problem for maritime safety and marine engineering.
Pioneering works by Munk et al. [1963]; Snodgrass et al. [1966] discovered a rich physics of the phenomenon and gave first examples of accurate measurements of magnitudes, periods and directional spreading of swell. Both articles [Munk et al., 1963;Snodgrass et al., 1966] contain thorough discussions of physical background of swell generation, attenuation and interaction with other types of ocean motions. Nonlinear wave-wave interactions have been sketched in these articles as a novelty introduced by the milestone papers by Phillips [1960] and Hasselmann [1962]. A possible important role of these interactions at high swell heights for relatively short time (fetch) of evolution has been outlined and estimated. The first estimates of the observed rates of swell attenuation have been carried out by Snodgrass et al. [1966] based on observation at near-shore stations. Their characteristic scale (e-folding) about 4000 km is consistent with some today results of the satellite tracking of swell [Ardhuin et al., 2009[Ardhuin et al., , 2010Jiang et al., 2016] and with treatment of these results within the model of swell attenuation due to coupling with turbulent atmospheric layer. Alternative model of turbulent wave flow attenuation [Babanin, 2006] predicts quite different algebraic law and stronger swell attenuation at shorter distances from the swell source [Young Copyright 2018 by the American Geophysical Union. 0148-0227/18/$9.00 et al., 2013]. It should be stressed that all the mentioned models treat swell as a quasi-monochromatic wave and, thus, ignore nonlinear interactions of the swell harmonics themselves and the swell coupling with locally generated windwave background. The latter effect can be essential as observations and simulations clearly show [Young, 2006;Badulin et al., 2008b, and refs. therein].
Today, linear models of swell are dominating. The swell is generally considered as a superposition of harmonics that do not interact to each other and, thus, can be described by the well-known methods of the linear theory of waves. Effects of weak dissipation, refraction by underwater relief and islands, sphericity of the Earth can be accounted for as well. Some features of the observed swell can be related to such models. For example, the observed linear growth of the swell frequency in a site can be explained as an effect of dispersion of a linear wave packet at long time. The linear models of swell dissipation lead, evidently, to exponential laws of the swell energy (height) attenuation and to attempts to estimate the e-folding scale from available experimental data [e.g. Snodgrass et al., 1966;Jiang et al., 2016].
Synthetic aperture radars (SAR) allow for spatial resolution up to tens meters and, thus, for detecting relatively long swell waves of a few hundred meters wavelength along their thousand miles tracks [e.g. Ardhuin et al., 2010;Young et al., 2013]. Satellite altimeters measure wave height averaged over a snapshot of a few square kilometers that is adequate to the today methods of statistical description of waves in the research and application models and also can be used for the swell tracking in combination with other tools [e.g. wave models as in Jiang et al., 2016]. Re-tracking of swell allows, first, for relating the swell events with their probable sources -stormy areas. Secondly, the swell transformation gives a clue to estimating effects of other motions of the ocean [e.g. Chen et al., 2002]. Such work requires adequate physical models of swell propagation and transformation as far as a number of parameters of sea environment remains beyond control.
The linear treatment of swell is not able to explain observed features of swell. First, the observed swell spectra exhibit frequency downshifting which is not predicted by linear or weakly nonlinear models of wave guide evolution [e.g. data of Snodgrass et al., 1966, and comments on these data by Henderson and Segur [2013]]. Secondly, these spectra show clearly invariance of their shaping that is unlikely to be appear in linear dispersive wave system. All the noted features are common for wave spectra described by the kinetic equation for water waves, the so-called Hasselmann [1962] equation.
In this paper we present results of extensive simulations of ocean swell within the Hasselmann equation. The simplest duration-limited setup has been chosen to obtain numerical solutions for the duration up to 10 6 seconds (11.5 days) for typical parameters of ocean swell (wavelengths 150 − 400 meters, initial significant heights 3 − 10 meters).
We analyze the simulation results from the viewpoint of the theory of weak turbulence [Zakharov et al., 1992]. The slowly evolving swell solutions appear to be quite close to the stationary milestone Kolmogorov-Zakharov solutions for water waves [Zakharov and Filonenko, 1966;Zakharov and Zaslavsky, 1982] in a frequency range and give estimates for the fundamental constants of the theory. We give a short theoretical introduction and present these estimates in the next section. In sect. 3 we relate results of simulations with properties of the self-similar solutions of the kinetic equation. Zaslavskii [2000] was the first who presented the self-similar solutions for swell assuming the angular narrowness of the swell spectra and gave a set of explicit analytical expressions for the swell. In fact, more general consideration in the spirit of Badulin et al. [2002Badulin et al. [ , 2005 leads to important findings and to a number of questions with no reference to the assumption of angular narrowness.
We demonstrate the well-known fact that is usually ignored: the power-law swell attenuation within the conservative kinetic equation. We show that it does not contradict to results of observations mentioned above. We also fix a remarkable feature of collapsing the swell spectra to an angular distribution that depends weakly on initial angular spreading. Such universality can be of great value for modelling swell and developing methods of its monitoring.
The paper is finalized by discussion of the effect wavewave interactions on swell attenuation. We consider a case where long-time swell evolution is perturbed by light (5 m/s) wind which direction is varying with inertial period. We show that the mechanism of inverse cascading due to weakly nonlinear wave-wave interactions provides an absorption of energy of locally generated short wind waves and an effective 'feeding' of long sea swell [Badulin et al., 2008b]. This effect remains beyond the attention of and cannot be treated within the models of swell attenuation mentioned above. Meanwhile, recent observations of swell from space show definitely possibility of swell amplification ['negative' dissipation in words of Jiang et al., 2016].

The Kolmogorov-Zakharov solutions
In this section we reproduce previously reported theoretical results on evolution of swell as a random field of weakly interacting wave harmonics. We follow the statistical theory of wind-driven seas [Zakharov , 1999] extending this approach to the sea swell which description within this approach is usually considered as questionable. A random wave field is described by the kinetic equation derived by Klauss Hasselmann in early sixties [Hasselmann, 1962] for deep water waves in absence of dissipation and external forcing Equation (1) is written for the spectral density of wave action N (k, x, t) = E(k, x, t)/ω(k) (E(k, x, t) is wave energy spectrum and wave frequency obeys linear dispersion relation ω = g|k|). Subscripts for ∇ corresponds to the twodimensional gradient operator in the corresponding space of coordinates x and wavevectors k (i.e. ∇x = (∂/∂x, ∂/∂y)).
The right-hand term S nl describes the effect of wave-wave resonant interactions and can be written in explicit form [see Badulin et al., 2005, for collection of formulas]. The cumbersome term S nl causes a number of problems for wave modelling where (1) is extensively used. Nevertheless, for deep water case one has a key property of the term homogeneity (2) that helps a lot in getting important analytical results. Stretching in κ times in wave scale or in ν times in wave action (κ, ν are positive) leads to simple re-scaling of the collision term S nl . This important property gives a clue for constructing power-law stationary solutions of the kinetic equation, i.e. solutions for the equation Two isotropic stationary solutions of (3) correspond to constant fluxes of wave energy and action in wave scales. The direct cascade solution [Zakharov and Filonenko, 1966] in terms of frequency spectrum of energy introduces the basic Kolmogorov constant Cp and describes the energy transfer to infinitely short waves with constant flux P . The wave action transfer to opposite direction of long waves is described by the inverse cascade solution [Zakharov and Zaslavsky, 1982] with wave action flux Q and another Kolmogorov's constant Cq.
An approximate weakly anisotropic Kolmogorov-Zakharov solution has been obtained by Katz and Kontorovich [1974] as an extension of (4) E (3) (ω, θ) = P 1/3 g 4/3 2ω 4 Cp + Cm gM ωP cos θ + . . . . (6) (6) associates the wave spectrum anisotropy with the constant spectral flux of wave momentum M and the so-called second Kolmogorov constant Cm. As it is seen from (6) the solution anisotropy vanishes as ω → ∞: wave spectra become isotropic for short waves. A harmony of the above set of the KZ solutions can be treated naturally within the dimensional approach (4-6): they are just particular cases of solutions of the form where G is a function of dimensional arguments scaled by spectral fluxes of energy P , wave action Q and wave momentum M . Originally, solutions (4-6) have been derived in quite sophisticated and cumbersome ways. Later on, simpler and more physically transparent approaches have been presented for generalization the set of these solutions and looking for their higher-order anisotropic extensions [Zakharov and Pushkarev , 1999;Balk , 2000;Pushkarev et al., 2003Pushkarev et al., , 2004Badulin et al., 2005;Zakharov , 2010]. In particular, the extension of the series in the right-hand side of (6) predicts the next term proportional to cos 2θ/ω 2 that is the second angular harmonics.
Swell solutions evolve slowly with time and, thus, give a good opportunity for discussing features of the KZ solutions (or, alternatively, the KZ solutions can be used as a benchmark for the swell studies). One of the key points of this discussion is the question of uniqueness of the swell solutions.
It can be treated in the context of general KZ solutions (7). While the principal terms of the general Kolomogorov-Zakharov solution corresponding to (4-6) have clear physical meaning of total fluxes of wave action (5), energy (4) and momentum (6) this is not the case of the higher-order terms. The link of these additional terms with inherent properties of the collision integral S nl or/and with specific initial conditions is not a trivial point. In this paper we give the very preliminary answer to this intriguing question.

Self-similar solutions of the kinetic equation
The homogeneity property (2) is extremely useful for studies of non-stationary (inhomogeneous) solutions of the kinetic equation. Approximate self-similar solutions for reference cases of duration-and fetch-limited development of wave field can be obtained under assumption of dominance of the wave-wave interaction term S nl [Pushkarev et al., 2003;Zakharov , 2005;Badulin et al., 2005;Zakharov and Badulin, 2011]. These solutions have forms of the so-called incomplete or the second type self-similarity [e.g. Barrenblatt, 1979]. In terms of frequency-angle dependencies of wave action spectra one has for the duration-and fetchlimited cases correspondingly [Badulin et al., 2005[Badulin et al., , 2007Zakharov et al., 2015] with dimensionless time τ and fetch χ Homogeneity properties (2) dictates 'magic relations' [in words of Zakharov , 2015, 2016] between dimensionless exponents pτ , qτ and pχ, qχ Dimensionless arguments of shaping functions Φp τ (ξ), Φp χ (ζ) contain free scaling parameters bτ , bχ Additional 'magic relations' coming from homogeneity property (2) fix a link between amplitude scales aτ , aχ and the bandwidth scales bτ , bχ of the self-similar solutions (8,9) Thus, 'magic relations' (11,13) reduce number of free parameters of the self-similar solutions (8,9) from four (two exponents and two coefficients) to two only: a dimensionless exponent pτ (pχ) and an amplitude of the solution aτ (aχ). The shaping functions Φ(ξ), Φ(ζ) in (8,9) require solution of a boundary problem for an integro-differential equation in self-similar variable ξ (or ζ for fetch-limited case) and angle θ [see sect. 5.2 Badulin et al., 2005, for details]. Simulations [e.g. Badulin et al., 2008a] show remarkable features of the shaping functions Φ(ξ), Φ(ζ). First, numerical solutions generally show relatively narrow angular distributions for Φp τ (ξ), Φp χ (ζ) with a single pronounced maximum near a spectral peak frequency ωp. It implies that the only one (or very few) of an infinite series of eigenfunctions of the boundary problem for the shaping functions Φp τ (ξ), Φp χ (ζ) contributes into wave spectra evolution in a wide range of initial and external forcing conditions. This treatment of the heavily nonlinear boundary problem in terms of a composition of eigenfunctions is possible in this case as demonstrated by Zakharov and Pushkarev [1999]. Two-lobe patterns are observed at higher frequencies (ω > 2ωp) in some cases as local maxima at oblique directions or as a 'shoulder' in wave frequency spectra. Their appearance is generally discussed as an effect of wind on wave generation [e.g. van Vledder , 2008, 2009].
Secondly, an important property of spectral shape invariance [terminology of Hasselmann et al., 1976] or the spectra quasi-universality [in words of Badulin et al., 2005] is widely discussed both for experimentally observed and simulated wave spectra. This invariance does not suppose a point-bypoint coincidence of properly normalized spectral shapes. Proximity of integrals of the shape functions Φp τ , Φp χ in a range of wave growth rates pτ , pχ appears to be sufficient for formulating a remarkable universal relationship for parameters of self-similar solutions (8,9) Here wave steepness µ is estimated from total wave energy E and spectral peak frequency ωp The 'number of waves' ν in a spatially homogeneous wind sea (i.e. for duration-limited wave growth) is defined as follows: For spatial (fetch-limited) wave growth the coefficient of proportionality C f in the equivalent expression ν = C f |kp|x (kp being the wavevector of the spectral peak) is close to the ratio between the phase and group velocities C ph /Cg = 2. A universal constant α0 ≈ 0.7 is a counterpart of the constants Cp, Cq of the stationary Kolmogorov-Zakharov solutions (4,5) and has a similar physical meaning of a ratio between wave energy and the energy spectral flux (in power 1/3). A remarkable feature of the universal wave growth law (14) is its independence on wind speed. This wind-free paradigm based on intrinsic scaling of wave development is capable to become a useful tool of analysis of both wind-wave growth and wind-free ocean swell .

Self-similarity of swell solutions
The self-similar solution for swell is just a member of a family of solutions (8,9) with particular values of temporal or spatial rates pτ = 1/11; qτ = 1/11 (17) pχ = 1/12; qχ = 1/12 (18) Exponents (17,18) provide conservation of the total wave action for its evolution in time (duration-limited setup) or in space (fetch-limited) On the contrary, total energy and wave momentum are formal constants of motion of the Hasselmann equation and decay with time t or fetch x The swell decay (20,21) reflects a basic feature of the kinetic equation for water waves: energy and momentum are not conserved [see Zakharov et al., 1992;Pushkarev et al., 2003, and refs. herein]. The wave action is the only true integral of the kinetic equation (1). The swell solution manifests another general feature of evolving wave spectra: the downshifting of the spectral peak frequency (or other characteristic frequency), i.e.
The universal law of wave evolution (14) is, evidently, valid for the self-similar swell solution as well with a minor difference in the value of the constant α0. As soon as this constant is expressed in terms of the integrals of the shape functions Φτ , Φχ and the swell spectrum shape differs essentially from ones of the growing wind seas this constant appears to be less than α0 of the growing wind seas. The theoretical facts presented above gives a background for analysis of results of our simulations.

Simulation setup
Simulations of ocean swell require special care. First of all, calculations for quite long time (up to 10 6 seconds in our case) should be accurate enough in order to fix relatively slow evolution of solutions and, thus, be able to relate results with the theoretical background presented above. We used the approach of our previous papers [Badulin et al., 2002[Badulin et al., , 2005[Badulin et al., , 2007[Badulin et al., , 2008a.
Duration-limited evolution has been simulated with the code based on WRT algorithm [Webb, 1978;Tracy and Resio, 1982] for different initial conditions (spatial spectra at t = 0). Frequency resolution for log-spaced grid has been set to (ωn+1 − ωn)/ωn = 1.03128266. It corresponds to 128 grid point in frequency range 0.02−1 Hz (approximately 1.5 to 3850 meters wave length). Thus, we used very fine frequency resolution in a wide range of wave scales [cf. Benoit and Gagnaire-Renou, 2007;Gagnaire-Renou et al., 2011] to trace rather slow evolution of swell. Standard angular resolution ∆θ = 10 • has been taken as adequate to goals of our study.
Initial conditions were similar in all series of simulations: wave action spectral density in a box of frequencies and angles was slightly (5% of magnitude) modulated in angles and had low (six orders less) background outside this box: (1 + 0.05 cos(θ 2 /2)), |θ| < Θ/2, ω l < ω < ω h 10 −6 N0, otherwise The modulations have been set in order to stimulate wavewave interactions as well as the collision integral S nl vanishes for N (k) = const. In contrast to wind waves where wind speed is an essential physical parameter that provides a useful physical scale, the swell evolution is determined by initial conditions only, i.e. by N0 and Θ in the setup (23). Dissipation was absent in the runs. Free boundary conditions were applied at the high-frequency end of the domain of calculations: generally, short-term oscillations of the Totally, more than 30 runs have been carried out for different initial conditions for duration 10 6 seconds. Below we focus on the series of Table 1 where initial wave heights were fixed (within 2%) and angular spreading varied from very narrow Θ = 30 • to almost isotropic Θ = 330 • (23). The frequency range of the initial perturbations was 0.1 − 0.4Hz.

Self-similarity features of swell
Evolution of swell spectra with time is shown in fig. 1 for the case sw330 of Table 1. The example shows a strong tendency to self-similar shaping of wave spectra. This remarkable feature has been demonstrated and discussed for    (15). The initial step-like spectrum evolves very quickly and keeps a characteristic shape for less than 1 hour. For 264 hours the spectral peak period reaches 11 seconds (the corresponding wavelength λ ≈ 186 meters) and wave steepness becomes µ = 0.023. The final significant wave height Hs ≈ 2.8 meters is essentially less than its initial value 4.8 meters. The solution parameters can be considered as typical ones for ocean swell.
Dependence of key wave parameters on time is shown in fig. 2 for different runs of the series of Table 1. Powerlaw dependencies of self-similar solutions (17,18,(20)(21)(22) are shown by dashed lines. In fig. 2a,b total wave energy E and the spectral peak frequency ωp show good correspondence to power laws of the self-similar solutions (8). All the run curves are approaching the theoretical dependencies in a similar way. By contrast, power-law decay of the x−component of wave momentum Mx depends essentially on angular spreading of initial wave spectra. While for narrow spreading (runs sw030 and sw050) there is no visible deviation from the law t −1/11 , wide-angle cases show these deviations definitely. The 'almost isotropic' solution for sw330 is tending quite slowly to the theoretical dependency in terms of wave momentum Mx (21). One can treat this transitional behavior by peculiarities of wave spectra relaxation from the 'almost isotropic' state to an inherent distribution with a pronounced anisotropy.
A simple quantitative estimate of the 'degree of anisotropy' is given in fig. 2d . Evolution of dimensionless parameter of anisotropy in terms of the approximate Kolmogorov-Zakharov solution (6) by Katz and Kontorovich [1974] is shown for all the cases of Table 1. We introduce parameter of anisotropy A as follows where total energy flux P (energy flux at ω → ∞) can be estimated from evolution of total energy  Table 1 (in  legend).
Similarly, total wave momentum provides an estimate of its flux in x-direction Spectral peak frequency ωp has been used for the definition of 'degree of anisotropy' A (24). Different scenarios are seen in fig. 2d depending on angular spreading of wave spectra. Nevertheless, a general tendency to a universal behavior at very large times (more than 10 6 seconds) looks quite plausible. Similar dispersion of different runs (different anisotropy of initial distributions) is seen in fig. 3 (14). One million seconds of the swell duration appear to be insufficient to present numerical arguments for the universality of the invariant (14) for the self-similar solutions (8). There is very likely a limit of the value at larger times. This limit is a bit less (by approximately 10%) than one for growing wind seas α0 ≈ 0.7. Again, the 'almost isotropic' solution shows its stronger departure from the rest of the series. The differences are better seen in angular distributions rather than in normalized spectral shapes ( fig. 4)  Spectrum at t = 0 is beyond of y-axis lower limit. function (or, more prudently, very few eigenmodes) survive in course of the swell evolution?
Normalized sections of spectra at the peak frequency ωp are shown in fig. 5 for the series runs at t = 10 6 seconds. 'The almost isotropic' run sw330 shows relatively high pedestal of about 3% of maximal value while other series have a background less than 0.2%. At the same time, the core of all distributions is quite close to a gaussian shape Evolution of angular spreading in time is shown in fig. 6 for 'the almost isotropic' run sw330 in absolute values. The sharpening of the angular distribution is accompanied by the peak increase of frequency-angle spectrum while the spectrum magnitude (averaged in angle) is slightly decaying with time (cf. fig. 1).
Angular spreading at higher frequency 1.25ω in fig. 7, first, illustrates a rapid saturation of high-frequency range: energy level at this frequency increases by 5 orders in magnitude for 1 hour only. Further evolution leads to pronounced sharpening of the distribution. The frequency 1.25ω is a characteristic one in the Discrete Interaction Approximation (DIA) for the collision integral S nl which is used extensively in research and forecasting models of sea waves [Hasselmann and Hasselmann, 1985]. Thus, fig. 7 likely warns possible problems in simulations of wave-wave interactions with the DIA: the numerical scheme should have rather fine distribution for resolving peculiar features of spectral functions. At large time this figure shows formation of side lobes: maxima of spectral densities for oblique counter-propagating harmonics. Such behavior of wave spectra is reported in many papers [e.g. Pushkarev et al., 2003;Bottema and van Vledder , 2008] and is usually discussed as an effect of wind input or as a transitional effect. Our simulations and the theoretical background of sect. 2.1 propose an alternative treatment  Figure 9. Left -estimate of the first Kolmogorov constant Cp, right -estimate of the second Kolmogorov constant Cm for the approximate anisotropic KZ solution (6). Time in hours is given in legend for the case sw170. of the effect in terms of properties of general stationary solutions (7) for the kinetic equation (1). Fig. 7 shows clearly the presence of the second angular harmonics that can be predicted within the formal procedure of Pushkarev et al. [2003Pushkarev et al. [ , 2004. Such treatment is not fully correct in the vicinity of the spectral peak but is still looks plausible. The second (and higher) harmonics of (7) is strongly decaying with frequency [e.g. Zakharov , 2010] and their accurate account requires special analysis. In the next section we present more evidences of their presence in the swell solutions.

Swell spectra vs KZ solutions
Very slow evolution of swell solutions in our simulations provides a chance to check relevance of the classic Kolmogorov-Zakharov solutions (4-7) to the problem under study. The key feature of the swell solution from the theoretical viewpoint is its 'hybrid' [in words of Badulin et al., 2005] nature: inverse cascading determines evolution of wave spectral peak and its downshifting while the direct cascading occurs at frequencies slightly (approximately 20%) above the peak. This hybrid nature is illustrated by fig. 8 for energy and wave momentum fluxes. In order to avoid ambiguity in treatment of the simulation results within the weak turbulence theory we will not detail this hybrid nature of swell solutions. Thus, general solution (7) in the form and its approximate explicit version (6) by Katz and Kontorovich [1971] will be used below for describing the direct cascading of energy and momentum at high (as compared to ωp) frequency. Positive fluxes P and M decays with time in good agreement with power-law dependencies (20) and have rather low variations in relatively wide frequency range 3ωp < ω < 7ωp in fig. 8. This domain of quasi-constant fluxes can be used for verification of results of the weak turbulence theory.
The first and the second Kolmogorov's constants can be easily estimated using the approximate solution (6) from combinations of along-and counter-propagating spectral densities as follows Cp = ω 4 (E(ω, 0) + E(ω, π)) 4g 4/3 P 1/3 (28) The corresponding values are shown in fig. 9. The found estimates of Kolmogorov's constants Cp ≈ 0.21 ± 0.01 and Cm ≈ 0.08 ± 0.02 are consistent with previous results [Lavrenov et al., 2002;Pushkarev et al., 2003;Badulin et al., 2005]. Note, that the cited papers used different definitions of Cp, Cm. Here we follow one of Zakharov [2010] (see eqs. 4-6). Estimate of the second Kolmogorov constant Cm requires special comments as far as its correctness can be affected by higher-order corrections (higher angular harmonics) of the solution (6). We see a possible effect of these corrections in relatively strong variations of estimates for the range ω/ωp < 4 in fig. 9b.
The consistency of the estimates of the Kolmogorov's constants for the swell solutions does not mean a good fit of the approximate weakly anisotropic solution (6) to the results of our simulations. At the same time, the robustness of the estimates provides a good benchmark for the balance of spectral fluxes and energy level in the wave spectra.

Discussion. Simulations for monitoring ocean swell 4.1. Swell attenuation within the kinetic equation
Results of the above consideration are two-fold. First, we illustrate the swell properties assuming the statistical approach for random water wave field (the Hasselmann equation) to be valid. We see that our simulations follow the basics of the weak turbulence theory quite well: deviations from self-similar asymptotics are generally small and have physical explanation. The second aspect comes from this fact of consistency of theory and simulations: we see an additional support of our approach for simulating the long  Figure 11. Results of simulations of swell under action of weak (5 m/s) rotating wind. Attenuation gives place to weak growth of swell after approximately 2 days of evolution. While dependence of wave height on time (a) is monotonous, wave momentum decay is accompanied by weak oscillations (b). Steps in peak wave period are from the simulation grid, no interpolation has been made. time swell evolution. Thus, the discussion of possible implications of our results to the problem of observation and monitoring ocean swell seems to be logical close of the work.
Dependence of wave height on time is shown in upper panel of fig. 10 for the 'almost isotropic' run sw330. Strong down up to 30% of initial value occurs at relatively short time of about one day. Essential part of the wave energy leakage corresponds to this transitional stage at the very beginning of swell evolution when swell is tending very rapidly to a self-similar asymptotics. Afterwards the decay becomes much slower following the power-law dependence for the selfsimilar solutions (20) fairly well.
For comparison with other models and available observations the duration-limited simulations have been recasted into dependencies on fetch through the simplest time-tofetch transformation [e.g. Hwang and Wang, 2004;Hwang, 2006]: The equivalent fetch is estimated as a distance covered by wave guide travelling with the group velocity of the spectral peak component. As seen in fig. 10 our model predicts an abrupt drop of wave height at distances shorter than 1000 km. In their turn, two quasi-linear models by Ardhuin et al. [2009] and Babanin [2006] predict gradual decay of swell up to 10% of initial wave height at distances 7000 km where our model shows qualitatively opposite weak attenuation. It should be noted that our model describes attenuation of the ocean swell 'on its own' due to wave-wave interactions without any external effects. Thus, the effect of an abrupt drop of wave amplitude at short time (fetch) should be taken into consideration above all others when discussing possible application of our results to swell observations and physical treatment of the experimental results.

Case study: swell under light wind
An accurate account of weak (as compared to the models of swell decay mentioned above) swell attenuation due to the effect of wave-wave interactions deserves a fuller explanation. Deficiency of the today understanding of swell physics and, in particular, swell attenuation can be illustrated by the results of our special case study. The setup of the case followed literally the one described in sect. 3.1 with a minor difference: the effect of light wind has been added. The wind input parameterization by Snyder et al. [1981] is used as a conventional one for the today wind-wave modelling. The experiment mimics an effect of wind that changes its direction with inertial period about 16 hours (inertial oscillations in the mid-latitude atmospheric boundary layer). The wind speed 5 m/s was approximately 3 times lower than phase speed of the spectral peak component at the start of the run (swell period Tp ≈ 10 seconds corresponds to phase speed higher than 15 m/s) and 5 times lower at the end of the run.
As seen from fig. 11 the swell decaying at initial stage started to grow after approximately 2 days of evolution when its period was about 15 seconds and wavelength exceeded 200 meters. While wave height Hs is growing ( fig. 11a) wave period Tp continue to follow the self-similar law t 1/11 (22) fairly well. All the dependencies in fig. 11 are monotonous except one of the wave momentum Mx (Tp in fig. 11d is not interpolated and is step-like because of numerical grid). The rotating wind is affecting the swell spectrum anisotropy but is feeding continuously wave energy. Wave-wave interactions are redistributing effectively energy of waves generated at high frequencies (3 or more times higher than peak frequency) and force swell to grow. This effect of the short wave absorption by swell has been discussed both for experimental data [Young, 2006] and results of simulations within the Hasselmann equation and weakly nonlinear dynamical equations for water waves [Badulin et al., 2008b]. Manifestations of the effect in global wave datasets [see Gulev and Grigorieva, 2003;Gulev et al., 2004, and http://www.sail.msk.ru/atlas/] have been analyzed and presented in terms of 'magic relations' (11) as well [Badulin and Grigorieva, 2012]. Recent attempts to re-track swell from satellite data and, thus, to estimate attenuation rates also show essential portion of 'negative dissipation rates' [in words of Jiang et al., 2016]. Most of the strange dissipation rates, in authors opinion, are 'not statistically significant'. At the same time, the distribution pattern of these rates correlates fairly well with results of analysis of the Voluntary Observing Ship (VOS) data [cf. figs 4,7 of Badulin and Grigorieva, 2012, and fig. 7 of Jiang et al. [2016]].
Thus, the presented case study shows deficiency of the today models of swell attenuation where the role of inherently nonlinear evolution of swell spectra is completely ignored. Swell is an effective 'devourer' of short wind waves [Badulin et al., 2008b] and it should be taken into account for the swell modelling.

Conclusions
We presented results of sea swell simulations within the framework of the kinetic equation for water waves (the Hasselmann equation) and treated these properties within the paradigm of the theory of weak turbulence. A series of numerical experiments (duration-limited setup, WRT algorithm) has been carried out in order to outline features of wave spectra in a range of scales usually associated with ocean swell, i.e. wavelengths larger than 100 meters and duration of propagation up to 10 6 seconds (≈ 11.5 days). It should be stressed that the exact collision integral S nl (nonlinear transfer term) has been used in all the calculations. Alternative, mostly operational approaches, like DIA (Discrete Approximation Approach) can corrupt the results qualitatively.
Fix the key results of the study 1. First, the classic Kolmogorov-Zakharov (KZ) isotropic and weakly anisotropic solutions for direct and inverse cascades are shown to be relevant to slowly evolving sea swell solutions. Estimates of the corresponding KZ constants are found to agree well with previous works. Thus, KZ solutions can be used as a benchmark for advanced approaches in the swell studies; 2. A strong tendency to self-similar asymptotics is observed. These asymptotics are shown to be insensitive to initial conditions. In particular, universal angular distributions of wave spectra at large times have been obtained for both narrow (appr. 30 • ) and almost isotropic initial spectra. The universality of the spectral shaping can be treated as an effect of mode selection when very few of eigenmodes of the boundary problem determines the system evolution. The inherent features of wave-wave interactions are responsible for this universality making effect of initial conditions insignificant. Initial conditions can affect rates of the tendency to this asymptotic universal state. Generally, we observe selfsimilar swell development in a background which is far from self-similar state; 3. We show that an inherent peculiarity of the Hasselmann equation, energy leakage, can also be considered as a mechanism of the sea swell attenuation. This mechanism is not accounted for by the today models of sea swell. In the meantime, the energy decay rates of sea swell in the numerical experiments, generally, do not contradict to results of swell observations. Moreover, ability of wave-wave interactions to redistribute energy and wave momentum over the whole range of wave scales can explain 'negative attenuation rates' found in satellite data.
The simulations of ocean swell within the durationlimited setup are just the very first step. Development of swell in space (fetch-limited setup) adds the effect of wave dispersion that can modify some results essentially. ----------------------