On ZRP wind input term consistency in Hasselmann equation

The new ZRP wind input source term (Zakharov et al. 2012) is checked for its consistency via numerical simulation of Hasselmann equation. The results are compared to field experimental data, collected at different sites around the world, and theoretical predictions of self-similarity analysis. Good agreement is obtained for limited fetch and time domain statements


Introduction
The scientific description of wind driven wave seas, inspired by solid state physics statistical ideas (see, for instance, Nordheim (1928)), was proposed by Hasselmann (1962Hasselmann ( , 1963 in the form of kinetic equation for waves: where ε = ε(ω k , θ, r, t) is the wave energy spectrum, as a function of wave dispersion ω k = ω(k), angle θ, two-dimensional real space coordinate r = (x, y) and time t. S nl , S in and S diss are nonlinear, wind input and wave-breaking dissipation source terms, respectively. Hereafter, only the deep water case, ω = √ gk is considered, where g is the gravity acceleration and k = | k| is the absolute value of wavenumber k = (k x , k y ).
Since Hasselmann's work, Eq.(1) has become the basis of operational wave forecasting models such as SWAN and Wavewatch III (Tolman 2013;SWAN), While the physical oceanography community consents on the general applicability of Eq.(1), there is no agreement on universal parameterizations of the source terms S nl , S in and S diss .
The S nl term was derived by different methods from the primodary Euler equations for free surface incompressible potential flow of a liquid by Hasselmann (1962Hasselmann ( , 1963 and Zakharov and Filonenko (1966). It is complex nonlinear operator acting on ε k , concealing hidden symmetries (Zakharov and Filonenko 1967;Zakharov et al. 1992). Resio and Perrie (1991) showed that those different forms are identical on the resonant surface ω k 1 + ω k 2 = ω k 3 + ω k 4 (2) The luxury of knowing the analytical expression for the S nl term, known in physical ocenography as XNL, is overshadowed by its computational complexity. Today, none of the operational wave forecasting models can afford performing XNL computations in real time.
Instead, an approximation, known as DIA and its derivatives, are used . The implication of such simplification is inclusion of a tuning coefficient in front of nonlinear term ; however, several publications have now shown that the form of the DIA does not provide a good approximation of the actual form of XNL . Consequently other source terms must be adjusted to allow the model Eq.
In contrast to S nl , the knowledge of the S in and S diss source terms is poor, so these now have many heuristic factors and terms included within them.
The creation of a reliable, well justified theory of S in has been hindered by strong turbulence presence in the air boundary layer over the sea surface. Even one of the most crucial elements of this theory, the vertical distribution of horizontal wind velocity in the region closest to the ocean surface, where wave motions strongly interact with atmospheric motions, is still the subject of debate. The history the development of different S in wind input forms is full of misconceptions, such as: the "fractional derivative method", wind driven waves studies in shallow basins, which fundamentally restrict the surface wave speed propagation, questionable data interpretations, and others. As a result, the values of different wind input terms scatter by the factor of 300 ÷ 500% (Badulin et al. 2005;Pushkarev and Zakharov 2016). Additional information on this detailed analysis of current state of wind input terms can be found in Pushkarev and Zakharov (2016).
Similar to the wind input term, there is little consent on the parameterization of the source dissipation term S diss . The physical dissipation mechanism, which most physical oceanographers agree on, is the effect of wave energy loss due to wave breaking, while there are also other dubious ad-hoc "long wave" dissipation source terms, having heuristically justified physical explanations. Wave breaking dissipation, known also as "white-capping dissipation", is an important physical phenomenon, not properly studied yet for the reasons of mathematical and technological complexity. Currently, there is not even agreement on the localization of wave breaking events in Fourier space. The approach currently utilized in operational wave forecasting models mostly relies on the dissipation, localized in the vicinity of the spectral energy peak. Recent numerical experiments show (Dyachenko et al. 2015;Zakharov et al. 2009), however, that such approach does not pass most of the tests associated with the essentially nonlinear nature of the Eq.(1).
Theoretical and experimental study of wave-breaking is extremely difficult task due to mathematical and technological complexities. Longuet-Higgins (1980b,a) achieved important results, but didn't accomplished the theory completely. Irisov and Voronovich (2011) studied wave-breaking of the short waves, "squeezed" by surface currents, caused by longer waves, and showed that they become steep and unstable. Our explanation is simpler, but has the same consequiences: the "wedge" formation, preceeding the wave breaking, causes the "fat tail" appearance is Fourier space. Subsequent smoothing of the tip of the wedge is equivalent to the "chopping off" of the developped high-frequency tail -the sort of natural low-pass filtering -leading to the loss of the wave energy. This effect is referred as "cigar cutting effect" (Pushkarev and Zakharov 2016). Both scenarios have the same consequences of wave surface smoothing, and are indirectly confirmed by presented numerical experiments.
Instead of following the previous path of time-consuming numerical and field experiments, the authors of the current manuscript are pursuing an alternative approach to definition of S in and S diss terms . Based on the leading nonlinearity role in Eq.(1) ( (Zakharov 2010;Zakharov and Badulin 2011)), it was desided to analyze a multi-parametric family of self-similar solutions of Eq.(1). The comparison with the results of field observations allowed to find a new wind input form, herein termed the ZRP S in wind input source term (Zakharov et al. 2012). The wave breaking dissipation term S diss has been chosen in the form of "implicit dissipation" via Phillips ∼ f −5 spectral continuation tail.
That framework reproduced the observations of a dozen of field experiments , i.e. selfsimilar exponents p = 1 and q = −0.3 of power dependencies of total wave energy ∼ χ p and spectral peak frequency ∼ ω q along the fetch (Pushkarev and Zakharov 2016).
In the following section we are describing the details of ZRP model and pertaining nu-merical results in time domain and duration limited statements.

Experimental evidence
Here we examine the empirical evidence from around the world, which has been utilized to quantify energy levels within the equilibrium spectral range by Resio et al. (2004). For convenience, we shall also use the same notation used by Resio et al. (2004) in their study, for the angular averaged spectral energy densities in frequency and wavenumber spaces: where f = ω 2π , α 4 is the constant, V is some characteristic velocity and β = 1 2 α 4 V g −1/2 . These notations are based on relation of spectral densities E(f ) and F (k) in frequency f = ω 2π and wave-number k bases: where c g = dω dk = 1 2·2π g f is the group velocity. The notations in Eqs.(4)-(5) are connected with the spectral energy density ǫ(ω, θ) through The Resio et al. (2004) analysis showed that experimental energy spectra F (k), estimated through averaging < k 5/2 F (k) >, can be approximated by linear regression line as a function of (u 2 λ c p ) 1/3 g −1/2 . Fig.1 shows that the regression line indeed, seems to be a reasonable approximation of these observations.
Here α 4 = 0.00553, u 0 = 1.93 m/sec, c p is the spectral peak phase speed and u λ is the wind speed at the elevation equal to a fixed fraction λ = 0.065 of the spectral peak wavelength 2π/k p , where k p is the spectral peak wave number. It is important to emphasize that Resio et al. (2004) experiments show that parameter β increases with development of the wind-driven sea, when f p decreases and C p increases. This observation is consistent with the weak turbulent theory, where β P 1/3 (Zakharov et al. 1992) (P is the wave energy flux toward small scales). Resio et al. (2004) assumed that the near surface boundary layer can be treated as neutral and thus follows a conventional logarithmic profile having Von Karman coefficient κ = 0.41, where z = λ · 2π/k p is the elevation equal to a fixed fraction λ = 0.065 of the spectral peak wavelength 2π/k p , where k p is the spectral peak wave number, and z 0 = α C u 2 ⋆ /g subject to Charnock (1955) surface roughness with α C = 0.015.

Theoretical considerations
Self-similar solutions of conservative kinetic equation were studied in Zakharov (2005), Badulin et al. (2005). In this chapter we study self-similar solutions of the forced kinetic equation where ǫ(ω, θ) = 2ω 4 g N(k, θ) is the energy spectrum. For our purposes, it is sufficient to simply use the dimensional estimate for S nl , Eq.(11) has a self-similar solution if where s is a constant. Looking for self-similar solution in the form we find The function F (ξ) has the maximum at ξ ∼ ξ p , thus the frequency of the spectral peak is The phase velocity at the spectral peak is According to experimental data, the main energy input into the spectrum occurs in the vicinity of the spectral peak, i.e. at ω ≃ ω p . For ω >> ω p , the spectrum is described by Here This integral converges if s < 2. For large ω More accurately Now supposing s = 4/3 and γ ≃ ω 7/3 , we get η = 1/3 , which is exactly experimental regression line prediction. Because it is known from regression line on Fig.1 that ξ = 1/3, we immediately get s = 4/3 and the wind input term At the end of the section, we present the summary of important relationships.
Wave action N, energy E and momentum M in frequency-angle presentation are: The self-similar relations for duration limited case: 9q − 2p = 1, p = 10/7, q = 3/7 s = 4/3 (29) The same sort of self-similar analysis gives self-similar relations for fetch limited case:

Numerical simulation
To check the self-similar conjecture Eq. (24), we performed the series of numerical simulations of Eq.(1) in the spatially homogeneous time domain ∂N ∂ r = 0 and spatially inhomogeneous fetch limited ∂N ∂t = 0 situations. The same wind input term Eq.(24) has been used in both cases in the form where u 10 is the wind speed, ρ air and ρ water are the air and water density correspondingly. It is conceivable to use more sophisticated expression for f (θ), for instance f (θ) = q(θ) − q(0).
The wind speed u 10 is taken here as the speed at a reference level of 10 meters. To make comparison with experimental results of Resio et al. (2004), we used relation u ⋆ ≃ u 10 /28 (see Golitsyn (2010)) in Eq.(9).
Both situations also need the knowledge of the dissipation term S diss , which is taken into account as was proposed by Resio et al. (2004), where white-capping was introduced implicitly through an f −5 energy spectral tail stretching in frequency range from f d = 1.1 to f max = 2.0. To date, this approach has been confirmed by both experimental observations and numerical experiments to provide an effective direct cascade energy sink at high frequencies.
The nonlinear S nl term has been calculated in the exact XNL form.

a. Time domain simulation
All time domain numerical simulations shown here have been started from uniform noise energy distribution in Fourier space. The check of consistency with the magic number (9q − 2p) (Pushkarev and Zakharov 2016), is presented on Fig.6. This also agrees with the self-similar prediction Eq.29. Fig.7 presents angle-integrated energy spectrum, as the function of frequency, in logarithmic coordinates. One can see that it consists of the segments of: • spectral maximum area • spectrum ω −4 • Phillips high frequency tail ω −5 The compensated spectrum F (k) · k 5/2 is presented on Fig.9. One can see plateau-like region responsible for k −5/2 behavior, equivalent to ω −4 tail in Fig.7. This exact solution of Eq.(10), known as KZ spectrum, was found by Zakharov and Filonenko (1967). One should note that most of the energy flux into the system comes in the vicinity of the spectral peak, as shown on Fig.8, providing significant intertial interval for KZ spectrum.
The angular spectral distribution of energy, presented on Fig.10, is consistent with the results of experimental observations that show a broadening of the angular spreading in both directions away from the spectral peak frequency.
Another important theoretical relationship, that can be derived from joint consideration of Eqs. (4), (6) and (22): presents theoretical equivalent of the experimental regression.
To compare the performance of our numerical simulation with the experimental analysis by Resio et al. (2004), presented in Fig.1, we show in Fig.11 plots of the compensated spectrum β = F (k) · k 5/2 as a function of (u 2 λ C p ) 1/3 /g 1/2 for wind speeds u 10 = 10.0 m/sec, along with the regression line from Resio et al. (2004). The overall correspondence is quite good.

b. Limited fetch numerical simulations
The limited fetch simulation was performed within the framework of the stationary version of the Eq.(1): where x is the coordinate orthogonal to the shore and θ is the angle between individual wavenumber k and the axis x.

Stationarity in Eq.(45) is somewhat difficult for numerical simulation, since it contains
the singularity in the form of cos θ in front of ∂ǫ ∂x . We overcame this problem of division by zero through zeroing one half of the Fourier space of the system for the waves propagating toward the shore. Since the energy in such waves is small with respect to waves propagating in the offshore direction, such approximation is quite reasonable for our purposes.
Since the wind forcing index s in the fetch-limited case is similar to that in the time domain, the numerical simulation of Eq.(45) has been performed for the same input functions as in the time doman case with the same low-level energy noise initial conditions in Fourier space.  The compensated spectrum F (k) · k 5/2 is presented on Fig.19. One can see plateu-like region responsible for k −5/2 behavior, equivalent to ω −4 tail in Fig.17. As in the time domain case, KZ solution (Zakharov and Filonenko 1967) also holds for limited domain case, and most of the energy flux into the system comes in the vicinity of the spectral peak as well, as shown on Fig.18, providing significant intertial interval for KZ spectrum. Angular spectrum distribution, consistent with the broadening of the spectrum in experimental observations is presented on Fig.20.   Fig.21 presents the plot of β = F (k) · k 5/2 as a function of (u 2 λ C p ) 1/3 /g 1/2 for wind speed u 10 = 10.0 m/sec, along with the regression line from Resio et al. (2004) and theoretical prediction Eq.(44). The overall correspondence is quite good.

Conclusion
We analysed the new ZRP form for wind input, proposed in Zakharov et al. (2012). Both numerical simulations for time domain and limited fetch cases, which used ZRP wind input term, XNL nonlinear term S nl and "implicit" high-frequency dissipation, shown good agreement with predicted self-similar properties of Hasselmann equation, experimentally obtained regression line Resio et al. (2004); Resio and Long (2007)and its theoretical prediction. The authors of the research hope that this new framework term will improve the quality of ocean wave prediction forecasts in operational models . Zakharov, V. E., A. O. Korotkevich, and A. O. Prokofiev, 2009: On dissipation function of ocean waves due to whitecapping. American Institute of Physics Conference Series, T. E. Simos, G.Psihoyios, and C. Tsitouras, Eds., Vol. 1168, 1229-1237. Zakharov, V. E., V. S. L'vov, and G. Falkovich, 1992: Kolmogorov Spectra of Turbulence I: Wave Turbulence. Springer-Verlag. Zakharov, V. E., D. Resio, and A. Pushkarev, 2012: New wind input term consistent with experimental, theoretical and numerical considerations. http://arxiv.org/abs/1212.1069/. Correlation of equilibrium range coefficient β with (u 2 λ c p ) 1/3 /g 1/2 based on data from six disparate sources. Adapted from Resio et al. (2004) 19  Resio et al. (2004); Resio and Long (2007). Crosses -results of numerical calculations for fetch limited case with wind speed U = 10 m/sec.