Balanced source terms for wave generation within the Hasselmann equation

The new Zakharov–Resio–Pushkarev (ZRP) wind input source term (Zakharov et al., 2012) is examined for its theoretical consistency via numerical simulation of the Hasselmann equation. The results are compared to field experimental data, collected at different sites around the world, and theoretical predictions based on self-similarity analysis. Consistent results are obtained for both limited fetch and duration limited 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 the Hasselmann equation (hereafter HE), also known as the kinetic equation for waves: where ε = ε(ω k , θ, r, t) is the wave energy spectrum as a function of wave dispersion ω k = ω(k), angle θ , twodimensional real space coordinate r = (x, y) and time t.S nl , S in and S diss are the 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 the vector wave number k = (k x , k y ).
Since Hasselmann's work, Eq. ( 1) has become the basis of operational wave forecasting models such as WAM, SWAN and Wavewatch III (Tolman, 2013;SWAN, 2015).While the physical oceanography community consents on the general applicability of Eq. ( 1), there is no consensus on universal parameterizations of the source terms S nl , S in and S diss .

The S nl term and weak turbulence theory
The HE is a specific example of the kinetic equation for quasi-particles, widely used in different areas of theoretical physics.There are standard methods for its derivation.In the considered case, two forms of the S nl term were derived by different methods from the Euler equations for free surface incompressible potential flow of a liquid by Hasselmann (1962Hasselmann ( , 1963) ) and Zakharov and Filonenko (1966).Resio and Perrie (1991) showed that they are identical on the resonant surface: (3) The S nl term is the complex nonlinear operator acting on ε k , concealing hidden symmetries (Zakharov and Filonenko, 1967;Zakharov et al., 1992) and cubic with respect to the spectrum ε.
Equation (1) sometimes is called "the Boltzmann equation" by the oceanographic community, although this is a misconception.The Boltzmann equation, derived in the nineteenth century for description of gas kinetics, is quadratic, rather than cubic, with respect to the distribution function.
Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.
V. Zakharov et al.: Balanced source terms for wave generation To understand the relation and difference between the Boltzmann equation and the HE, one should recall the abovementioned Nordheim (1928) equation.This equation, applicable to quantum quasi-particles, contains both quadratic and cubic terms.Hence, the Boltzmann equation and the HE present opposite limiting cases of a general quantum kinetic equation.
Purely cubical (applicable to classical waves, not to classical particles) systems are relatively new objects in physics.Such equations describe the simplest case of the wave turbulence by the theory called weak turbulence theory (WTT) (Zakharov et al., 1992).
It is unfortunate that the discussion of the HE in the context of WTT has been overlooked by a major part of the oceanographic community for many years now.The community accepts, nevertheless, the HE as the basis for the operational wave forecasting models, thereby believing de facto in WTT without fully appreciating its ramifications.
The WTT essentially differs from the kinetic theory of classical particles and quantum quasi-particles.In the "traditional" gas kinetics (both classical and quantum) the basic solutions are thermodynamic equilibrium spectra, such as Boltzmann and Planck distributions.In the WTT such solutions, though formally existing, play no role -they are non-physical.The physically essential solutions are the nonequilibrium Kolmogorov-Zakharov spectra (or KZ spectra, Zakharov et al., 1992), which are the solutions of the corresponding kinetic equation The simplest one is the Zakharov-Filonenko (hereafter ZF) solution (Zakharov and Filonenko, 1966), which is the subclass of KZ solutions: where P is the energy flux toward high wave numbers.
The accuracy advantage of knowing the analytical expression for the S nl term, also known in physical oceanography as the XNL, is overshadowed by its computational complexity.Today, none of the operational wave forecasting models can afford to perform XNL computations in real time.Instead, the operational approximation, known as DIA and its derivatives, is used to replace this source term.The implication of such simplification is the inclusion of a tuning coefficient in front of the nonlinear term; however, several publications have shown that DIA does not provide a good approximation of the actual XNL form.The paradigm of replacement of the XNL by DIA and its variations leads to even more grave consequences: other source terms must be adjusted to allow the model Eq.(1) to produce desirable results.In other words, deformations suffered by the XNL model due to the replacement of S nl by its surrogates need to be compensated for by non-physical modification of other source terms to achieve reasonable model behavior in any specific case, leading to a loss of physical universality in the HE model.

Operational formulations for the wind energy input
S in and wave energy dissipation S diss terms In contrast to S nl , the knowledge of S in and S diss source terms is poor; furthermore, both include many heuristic factors and coefficients.The creation of a reliable, well-justified theory of S in has been hindered by strong turbulent fluctuations, uncorrelated with the wave motions, in the 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 of the development of different wind input forms is full of heuristic assumptions, which fundamentally restrict the magnitude and directional distribution of this term.As a result, the values of different wind input terms scatter by a factor of 300-500 % (Badulin et al., 2005;Pushkarev and Zakharov, 2016).For example, experimental determination of S in , as provided by direct measurements of the momentum flux from the air to the water, cannot be rigorously performed in a laboratory due to gravity wave dispersion dependence on the water depth, as well as problems with scale effects in laboratory winds.Additional information on the detailed analysis of the current state of the art of wind input terms can be found in Pushkarev and Zakharov (2016).Similar to the wind input term, there is little consensus 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 with heuristically justified physical explanations.Currently, there is not even an agreement on the location 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 (Pushkarev and Zakharov, 2016;Dyachenko et al., 2015;Zakharov et al., 2009), however, that such an approach does not pass most of the tests associated with the essentially nonlinear nature of the HE Eq. (1).The next chapters present a balanced set of wind energy input and wave energy dissipation source terms, based on WTT and experimental data analysis.Further, they are numerically checked to comply with WTT predictions and experimental observations.As mentioned above, contrary to previous attempts to build the detailed-balance source terms, the current approach is based neither on the development of a rigorous analytic theory of turbulent atmospheric boundary layers nor on reliable and repeatable air to ocean momentum measurements.The new S in is constructed in the artificial way realizing, in a sense, "the poor man approach", based on the finding of a two-parameter family of HE self-similar solutions and their restriction to the single-parameter one with the help of comparison with the data of experimental observations, accumulated for several decades.Section 2 presents experimental evidence of wave energy spectrum characteristics in the form of a specific regression line, found by Resio et al. (2004).The analytic form of this regression line will play a crucial role in narrowing the circle of possible outcomes, obtained with WTT analysis.
Section 3 studies self-similar solutions of the HE -a kinetic equation for surface ocean waves, starting with the analysis of the behavior of the dissipationless HE in infinite space, containing the wind source term in power function form.This approach is similar in spirit to one realized by Zakharov and Filonenko (1967) for finding the solution of the equation S nl = 0 in the infinite Fourier domain, which derived the ZF spectrum ε ∼ P 1/3 ω 4 , where P is the energy flux toward high wave numbers.The Fourier domain in both situations does not contain any dissipation function: its role is played by infinite phase volume as the effective energy sink at infinitely high wave numbers.
Such a situation is similar to one realized in incompressible liquid turbulence for large Reynolds numbers, where the energy distribution is given by the famous Kolmogorov spectrum, transferring the energy from large to small scales, where the energy dissipation is realized due to viscosity, but the viscosity coefficients, i.e., the dissipation details, are not included in the final Kolmogorov spectrum expression.The ZF spectrum and its further KZ generalizations are in this sense the ideological Kolmogorov spectrum counterparts, having the significant difference that the Kolmogorov spectrum is a plausible conjecture, while KZ spectra are the exact solutions of the wave kinetic equation.
Since the current research is application oriented, it is important to understand why this formally academic approach is connected with reality.In this context, there is no such thing as the dissipation at infinitely small waves in nature: however, it is clear that the existence of an absorption at sufficiently high finite frequencies provides a wave scale in real applications that still preserves the KZ solutions, found from the HE equation in infinite space.
As was mentioned before, this statement was confirmed in a different physical context with radically different inertial ranges (the wave-number band between characteristic energy input and characteristic wave energy dissipation), showing KZ solutions with different corresponding indices.As for the considered case of gravity waves on the surface of a deep fluid, KZ spectra have been routinely observed in multiple experiments.The results, published before 1985, are summarized by Phillips (1985).Thereafter, they were observed and discussed by Long and Resio (2007).A complete survey of all measurements requires a separate comprehensive paper, which is in our plans for the future.
The assumed close relation of the HE in the infinite space and finite domain, bounded by high-frequency dissipation, also has a much deeper meaning, consisting in the fact that S nl is the leading term of the HE (Zakharov, 2010;Zakharov and Badulin, 2011).This allows further use of the solutions found from the "zero-dissipation" HE Eq. ( 13) in infinite space for "practical" Fourier domains with the dissipation localized at finite high enough wave numbers.They take the form of a two-parameter family of self-similar solutions, which can be further restricted to the single-parameter one using experimental regression dependence, presented in Sect. 2. These self-similar solutions present realistic HE solutions and describe a broad class of wave-energy spectra observed in ocean and wave-tank experiments.
The indices, corresponding to self-similar solutions, allow one to wrap up Sect. 3 with the specific form of the wind input term in infinite phase space, called the ZRP wind input term (Zakharov et al., 2012;Pushkarev and Zakharov, 2016), with an arbitrary coefficient in front of it.Now, the theoretical part of the wind source term S in construction is finished, but the obtained model is not suitable yet for numerical simulation, since to perform in finite phase space, it has to be augmented with the wave-breaking dissipation term.
Section 4 explains the dissipation function used in the presented model.The wave-breaking dissipation, also known as "white-capping dissipation", is an important physical phenomenon not properly studied yet for reasons of mathematical and technological complexity.Longuet-Higgins (1980a, b) achieved important results, but did not accomplish the theory completely.Irisov and Voronovich (2011) studied the wave-breaking of 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 consequences: the "wedge" formation, preceding the wave breaking, causes the "fat tail" appearance in Fourier space.Subsequent smoothing of the tip of the wedge is equivalent to a "chopping off" of the developed high-frequency tail in Fourier space -a sort of natural low-pass filtering -leading to the loss of the wave energy.Both scenarios lead to smoothing of the wave surface, and are indirectly confirmed by the numerical experiments presented in the current study.

V. Zakharov et al.: Balanced source terms for wave generation
There is considerable freedom in choosing a specific analytic form of such a high-frequency dissipation term, given the lack of a generally accepted rigorous derivation for this mechanism.Consequently, one can choose a preferred one and possibly justify it, but any particular choice will be questioned since it will remain somewhat artificial.Because of that, our motivation was that at the current stage of development, we considered simplicity as a primary motivating factor.Instead of following the previous path of time-consuming numerical and empirical formulations based on field experiments, the authors decided to continue the spectrum from some specific frequency point, well above the spectral peak, with the Phillips law ∼ ω −5 , which decays faster than the equilibrium spectrum ω −4 and therefore corresponds to a net wave energy absorption.Although a version of this concept was incorporated by Janssen (2009), detailed forms of this source term have not been developed to date, other than that the spectrum at high frequencies appears to consistently tend toward an ∼ ω −5 form as noted by Phillips (1985).Additional evidence for a transition from ∼ ω −4 to ∼ ω −5 at frequencies above the equilibrium range comes from analysis of multiple data sets by Resio et al. (2004).In that paper the transition from ∼ ω −4 to ∼ ω −5 occurs approximately at f d = 1.1 Hz; i.e., the physical spectrum has to be continued from this point by ∼ ω −5 .
The spectrum amplitude at the junction frequency f d is dynamically changing in time.It is important that this analytic continuation contributes to a differential in inverse action, which also affects frequencies lower than f d , since the nonlinear interaction term S nl is calculated over both "dynamic" and fixed Phillips areas.Therefore, the Phillips part of the spectrum "sends" the information about the presence of the dissipation above f d to the rest of the spectrum.
At this point, all that remains for source-term closure in the HE model is the coefficient in front of the wind input term, since it is not well defined experimentally.If we carry out the numerical simulation with some arbitrary chosen coefficient, we could obtain a range of spectral energies but would retain the qualitative properties of the HE, like the ∼ ω −4 spectrum, spectral peak down-shift and peak frequency behavior in accordance with self-similar laws.
To solve this, we choose the wind source coefficient to reproduce the same wave energy growth as was observed in field experiments.The value of this coefficient, found from the comparison with field observations of wave energy growth, is equal to 0.05.This step completes the construction of the HE model.
In the next sections we proceed with numerical simulations based on the HE model described above.Section 5 discusses the details of the numerical model setup.Section 6 describes the duration limited numerical simulation, which is the subject of more academic than applied interest, targeted at self-similarity concept support, while the limited fetch numerical simulation results, described in Sect.7, besides academic interest, are the subject of comparison with the field experiments.A check of the compliance of numerical results with field experimental measurements is presented in Sect.8.

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 g f is the group velocity.The notations in Eqs. ( 6) and ( 7) 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 a linear regression line as the function of (u 2 λ c p ) 1/3 g −1/2 .Figure 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 s −1 , 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 the 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); here 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

Theoretical considerations
Self-similar solutions consistent with the conservative kinetic equation were studied in Zakharov (2005) and Badulin et al. (2005).
In this section we study self-similar solutions of the forced kinetic equation where (ω, θ ) = 2ω 4 g N (k, θ ) is the energy spectrum.One should note that this equation does not contain any explicit wave dissipation term; the role of dissipation is played by the existence of the energy sink at infinitely high wave numbers, in the spirit of the WTT; see Zakharov and Filonenko (1967) and Zakharov et al. (1992).
For our purposes, it is sufficient to simply use the dimensional estimate for S nl , Eq. ( 13) has a self-similar solution if where s is a constant.Looking for a self-similar solution in the form we find The function F (ξ ) has a 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 the Zakharov-Filonenko tail 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 the experimental regression line prediction.
Because it is known from the regression line in Fig. 1 that ξ = 1/3, we immediately get s = 4/3 and the wind input term For many years, the assumption has been that there could be a net input or dissipation within the equilibrium range;  Resio et al. (2004) suggest that the existence of significant net energy input or dissipation within the frequency range would tend to force the spectrum away from an f −4 form, contrary to the pattern found in field measurements.If we assume that the wind source is primarily centered on the spectral peak, the only missing component in our numerical solution is an unknown coefficient in front of it, which will be defined later from the comparison with total energy growth in experimental observations.Another important theoretical relationship that can be derived from joint consideration of Eqs. ( 6), ( 8) and ( 24) is which shows a theoretical equivalence to the experimental regression, where λ is an unknown constant, defined experimentally.
At the end of the section, we present the summary of important relationships.
Wave action N , energy E and momentum M in frequencyangle presentation are The self-similar relations for the duration limited case are given by = t p+q F (ωt q ), (31) The same sort of self-similar analysis gives self-similar relations for the fetch limited case: 4 The details of "implicit" dissipation Now that the construction of the ZRP wind input term with the unknown coefficient has been accomplished in the spirit of WTT in the previous chapter, the HE model, suitable for numerical simulation, still misses the dissipation term localized at finite wave numbers -there is no such thing as the infinite phase volume in reality: the real ocean Fourier space is confined by a characteristic wave number corresponding to the start of the dissipation effects caused by the wavebreaking events.
There is a lot of freedom in choosing the dissipation term.Since there is no current interpretation of the wave-breaking dissipation mechanism, one can choose it in whatever shape it is preferred, but any particular choice will be questioned since it is an artificial one.
Because of that, the motivation consisted in the fact that at the current "proof of concept" stage one needs to know the effective sink with the simplest structure.Continuation of the spectrum from ω d with the Phillips law A(ω d ) • ω −5 (see Phillips, 1966), decaying faster than the equilibrium spectrum ω −4 , will get high-frequency dissipation.The corresponding analytic parameterization of this dissipation term will be unknown, while not in principle impossible to figure out in some way.One should note that this method of dissipation is not our invention: it is described in Janssen (2009).
Specifically, the coefficient A(ω d ) in front of ω −5 is unknown but is not required to be defined in an explicit form.Instead, it is dynamically determined from the continuity condition of the spectrum, at frequency ω d , on every time step.In other words, the starting point of the Phillips spectrum coincides with the last point of the dynamically changing spectrum, at the frequency point ω d = 2π f d , where f d 1.1 Hz, as per Long and Resio (2007).This is the way the high-frequency "implicit" damping is incorporated into the alternative computational framework of the HE.The question of the finer details of the high-frequency "implicit" damping structure is of secondary importance, at the current "proof of concept" stage.
The whole set of the input and dissipation terms is accomplished now with one uncertainty: the explained approach leaves one parameter arbitrary -the constant in front of the wind input term.We choose it to be equal to 0.05 from the condition of the reproduction of the field observations of wave energy growth along the fetches, analyzed in Badulin et al. (2007).

Numerical validation of the relationship
To check the self-similar hypothesis posed in Eq. ( 26), we performed a series of numerical simulations of Eq. (1) in the spatially homogeneous duration limited ∂N ∂r = 0 and spatially inhomogeneous fetch limited ∂N ∂t = 0 situations.All simulations used the WRT (Webb-Resio-Tracy) method (see Tracy and Resio, 1982), which calculates the nonlinear interaction term in the exact form.The presented numerical simulation utilized the version of the WRT method, previously used in Webb (1978), Resio and Perrie (1989), Perrie and Zakharov (1999), Pushkarev et al. (2003), Long and Resio (2007), Korotkevich et al. (2008), Zakharov and Badulin (2015), and Pushkarev and Zakharov (2016), and used the grid of 71 logarithmically spaced points in the frequency range from 0.1 to 2.0 Hz and 36 equidistant points in the angle domain.The constant time step in the range between 1 and 2 s has been used for explicit first-order accuracy integration in time.
There is a balance between the number of nodes of the grid and the volume of the calculation to be performed.The particular version of the WRT model has been tuned to the minimum grid number of nodes to solve realistic physical problems, but is still fast enough to simulate them over a reasonable time span.The correctness of this statement is confirmed by the multiple numerical experiments cited above, reproducing mathematical properties of the Hasselmann equation.
For convenience, we present the pseudo-code used for the main cycle of the described model.

Return to step 1.
All numerical simulations discussed in the current paper have been started from a uniform noise energy distribution in Fourier space ε(ω, θ ) = 10 −6 m 4 , corresponding to a small initial wave height with an effectively negligible nonlinearity level.The constant wind of speed 10 m s −1 was assumed to blow away from the shoreline, along the fetch.The assumption of constant wind speed is a necessary simplification due to the fact that the numerical simulation is being compared to various data from field experiments, and the considered setup is the simplest physical situation which can be modeled.
The same ZRP wind input term Eq. ( 26) has been used in both cases as where U is the wind speed at the reference level of 10 m, and ρ air and ρ water are the air and water density, respectively.It is conceivable to use a more sophisticated expression for q(θ ), for instance q(θ ) = q(θ ) − q(0).To make direct comparison with experimental results of Resio et al. (2004), we used the relation u * U/28 (see Golitsyn, 2010) in Eq. ( 11).Frequencies f min and f d depend on the wind speed and should be found empirically.In current numerical experiments for U = 10 and U = 5 m s −1 , f min = 0.1 Hz and f d = 1.1 Hz.This choice is justified by the obtained numerical results.The above-described "implicit dissipation" term S diss has played the dual role of a direct energy cascade flux sink due to wave breaking as well as a numerical scheme stabilization factor at high wave numbers.

Duration limited numerical simulation
The duration limited simulation has been performed for a wind speed of U = 10 m s −1 .
Figure 2 shows the total energy growth as the function of time, consistent with the self-similar prediction Eq. ( 34) for index p = 10/7, supplied with the empirical coefficient in front of it; see Fig.One should specifically elaborate on the local index p numerical calculation procedure for Fig. 3. First, the total energy function was smoothed via a moving average, then the corresponding derivative is estimated numerically via finite differences, and finally a moving average is used to obtain the time-varying index value.32) is represented by the horizontal dashed line.
The relatively small systematic deviation from self-similar behavior, visible in Fig. 2, is connected with the following two facts.
First, the transition process in the beginning of the simulation, when the wave system behavior is far from a self-similar one.The self-similar solution is a pure power function, which does not take into account the initial transition process, and which causes the systematic difference.This systematic difference could be diminished via a parallel shift, which would take into account the initial transition process.Such a parallel shift is equivalent to starting the simulation from different initial conditions.
The second fact is the asymptotic nature of the self-similar solution, producing an evolution of the simulated wave system toward self-similar behavior with increasing time.As seen in Fig. 3, the numerical value of the local exponent converges to the theoretical value p = 10/7, reaching approximately 6 % accuracy for a sufficient dimensionless time 3 × 10 4 .
The dependence of the mean frequency on time, shown in Fig. 4, is consistent with the self-similar dependence found in Eq. ( 36) for q = 3/7, supplied with the empirical coefficient in front of it: see Fig. 5.
The systematic deviation of two lines in Fig. 4 remains within 3 % of the target value q = 3/7 for the same reasons as for wave-energy behavior -the transition process in the beginning of the simulation and the asymptotic nature of the self-similar solution.
A check of the consistency with the "magic number" 9q − 2p = 1 (see Eq. 32) is presented in Fig. 6.The reason for systematic deviation from the target value 1 is obviously connected with the reasons for the systematic deviations of p and q, as the "magic number" is calculated as their linear combination, reaching the accuracy of approximately 10 % for a long enough dimensionless time of 3 × 10 4 .
One should note that indices p and q and the "magic relation" 9q − 2p exhibit asymptotic convergence to the corresponding target values.λ c p ) 1/3 /g 1/2 .Dashed line -theoretical prediction Eq. ( 10) for λ = 2.74; dotted line -experimental regression line from Resio et al. (2004) and Long and Resio (2007).Line connected diamonds -results of numerical calculations for the wind speed U = 10 m s −1 duration limited case.Being parameterized by dimensionless time tg/U , the numerical simulation trajectory evolves from the left to the right on the graph, covering a time span from tg/U = 0 to tg/U 3.5 × 10 5 .
One can see a plateau-like region responsible for k −5/2 behavior, equivalent to the ∼ f −4 tail in Fig. 7.This shape of the spectrum is similar to that observed by Long and Resio (2007).This exact solution of Eq. ( 12), known as the KZ spectrum, was found by Zakharov and Filonenko (1967).The universality of f −4 asymptotic for the "inertial" (also known as "equilibrium" in oceanography) range between spectral peak energy input and high-frequency energy dissipation areas has been observed in multiple experimental field observations and is accepted by the oceanographic community after the seminal work of Phillips (1985).One should note that most of the energy flux into the system comes in the vicinity of the spectral peak, as shown in Fig. 9, providing a significant inertial interval for the KZ spectrum.
The angular spectral distribution of energy, presented in Fig. 10, is consistent with the results of experimental observations (Resio et al., 2011) that show a broadening of the angular spreading in both directions away from the spectral peak frequency.
To compare the duration limited numerical simulation results with the experimental analysis by Resio et al. (2004), presented in Fig. 1, Fig. 11 shows the function β = F (k) • k 5/2 as the function of (u 2 λ C p ) 1/3 /g 1/2 for wind speed U = 10 m s −1 , along with the regression line from Resio et al. (2004) and the theoretical prediction Eq. ( 27) for λ = 2.74.The numerical results and theoretical prediction line fall within a very small rms deviation (r 2 = 0.939; see Fig. 1) from the regression line.One should note asymptotic convergence of the numerical simulation results to the theoretical line.

Limited fetch numerical simulation
The limited fetch simulation was performed in the framework of the stationary version of Eq. ( 1): where x is chosen as the coordinate axis orthogonal to the shore and θ is the angle between the individual wave number k and the axis x.To find the dependence on the wind speed, directed off the shore, two numerical simulations for wind speeds of U = 5 m s −1 and U = 10 m s −1 have been performed.
The stationarity in Eq. ( 47) is somewhat difficult for numerical simulation, since it contains a singularity in the form of cos θ in front of ∂ ∂x .This problem was overcome by 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 an approximation is quite reasonable for our purposes.
Since the wind forcing index s in the fetch limited case is similar to that in the duration limited case, the numerical simulation of Eq. ( 47) has been performed for the same input functions as in the duration limited case with the same lowlevel energy noise initial conditions in Fourier space.
Figure 12 presents total energy growth as a function of fetch, consistent with the self-similar solution Eq. ( 40 index p = 1, and its appropriate empirical coefficient.The corresponding values of indices p along the fetch are presented in Fig. 13.The small amplitude oscillations observed in the index behavior can be attributed to the finite grid resolution used in the simulation. The wave evolution for the wind speed U = 5 m s −1 case is expected to be slower than for the U = 10 m s −1 case due to the weaker nonlinear interaction term.One can see, indeed, slower asymptotic convergence of the calculated total , where n(ω, θ) = ε(ω,θ) ω is the wave action spectrum, for wind speed 10 m s −1 (solid line) and 5 m s −1 (dashed line).The dash-dotted line is the self-similar dependence 3.4 • ( xg U 2 ) −0.3 with the empirical coefficient in front of it.
energy local power index to the target value p = 1 for the U = 5 m s −1 case compared to the U = 10 m s −1 case.The deviation of results from the U = 10 m s −1 case relative to the target value does not exceed an error of about 5 %, while for the U = 5 m s −1 case the error does not exceed 20 %.The role of the relatively short (in time) non-self-similar development of the wave system at the very beginning of the fetch should be noted as well as the factor contributing to the deviation from the target value of index p = 1: the wave system obviously needs some time to evolve into a fully self-similar mode.
The dependence of the mean frequency on the fetch, shown in Fig. 14, is consistent with the self-similar dependence Eq. ( 42) for index q = 0.3, supplied with the empirical coefficient in front of it.The small amplitude oscillations observed in index behavior can be attributed to the finite grid resolution used in the simulation, since the spectral peak moves continuously between discrete frequencies in a manner that cannot be matched in these discretized simulations.
The local values of indices q for two different wind speed amplitudes are presented in Fig. 15 along with the target value of the self-similar index q = 0.3.After sufficient fetch one can see only about 14 % deviation from the target value for the U = 10 m s −1 case and about 2.5 % for the U = 5 m s −1 case.
The reasons for the 10 % systematic deviation from the self-similar solutions observed in the lines in Fig. 14 behavior -the transition process in the beginning of the simulation and the asymptotic nature of the self-similar solution.
The check on the consistency of the calculated "magic number" (10q − 2p) (see Eq. 38) is presented in Fig. 16.The reason for systematic deviation from the target value 1 is obviously connected with the systematic deviations of p and q, as the "magic number" is calculated as their linear combination, reaching the accuracy of approximately 10 % for fetches 3×10 4 .As noted previously, the small amplitude oscillations observed in the indices' behavior can be attributed to the finite grid resolution used in the simulation.
One should note that indices p and q and the "magic relation" 10q − 2p exhibit asymptotic convergence to the corresponding target values.
Figure 17 presents an angle-integrated energy spectrum, as the function of frequency, in logarithmic coordinates.As  λ c p ) 1/3 /g 1/2 .Dashed line -theoretical prediction Eq. ( 10) for λ = 2.11; dotted line -experimental regression line from Resio et al. (2004) and Long and Resio (2007).Line connected diamonds -the results of numerical calculations for the wind speed U = 10 m s −1 limited fetch case.Being parameterized by the dimensionless fetch coordinate χ = xg U 2 , the numerical simulation trajectory evolves from the left to the right on the graph, covering a fetch span from χ = 0. to χ 3.0 × 10 4 .k −5/2 behavior, equivalent to the ω −4 tail in Fig. 17 and similar to that observed by Long and Resio (2007).As in the duration limited case, the KZ solution (Zakharov and Filonenko, 1967) also holds for the fetch limited case, and most of the energy flux into the system comes in the vicinity of the spectral peak as well, as shown in Fig. 18, providing a significant inertial (equilibrium) range for the KZ spectrum  Young (1999).The original caption is "A composite of data from variety of studies showing the development of the non-dimensional energy, ε as a function of non-dimensional fetch, χ .The original JONSWAP study (Hasselmann et al., 1973) used the data marked, JONSWAP, together with that of Burling (1959) andMitsuyasu (1968).Also shown are a number of growth curves obtained from the various data sets.These include JONSWAP Eq. (5.27), Donelan et al. (1985) Eq. (5.33) and Dobson et al. (1989) Eq. (5.38)." between spectral peak energy input and high-frequency energy dissipation areas.
The angular spectral distribution of energy, presented in Fig. 21, as in the duration limited case, is consistent with the results of experimental observations by Resio et al. (2011) that show a broadening of the angular spreading in both directions away from the spectral peak frequency.
The excess spectral energy at very oblique angles is a numerical artifact connected with the specifics of how the limited fetch statement is simulated here, i.e., the abovementioned singularity presence on the left-hand side of Eq. ( 47) at θ = ±π/2.
The detailed structure of angular spreading for both the duration and limited fetch cases is given in Fig. 19.The time that would be required to produce such a pattern is far in excess of the time for this excess energy to be removed from the equilibrium range by the nonlinear flux and can be shown to vanish when a time-space simulation is used instead of the stationary solution assumed here.
It is clearly seen that the "blobs" in the limited fetch case contain no more than 5 % of the total energy of the corresponding spectrum and could be neglected for the purposes of the presented research.
To compare the limited fetch numerical simulation results with the experimental analysis by Resio et al. (2004), presented in Fig. 1, Fig. 22 shows the function β = F (k)•k 5/2 as a function of (u 2 λ C p ) 1/3 /g 1/2 for wind speed U = 10.0 m s −1 , along with the regression line from Resio et al. (2004) and its theoretical prediction Eq. ( 27) for λ = 2.11.The numerical results and theoretical prediction line fall within the rms deviation (r 2 = 0.939; see Fig. 1) from the regression line.One should note asymptotic convergence of the numerical simu-  Young (1999).The original caption is "A composite of data from a variety of studies showing the development of the non-dimensional peak frequency, ν as a function of non-dimensional fetch, χ .The original JONSWAP study (Hasselmann et al., 1973) used all the data shown with the exception of that marked Donelan et al. (1985) and Dobson et al. (1989).Also shown are a number of growth curves obtained from the various data sets.These include JONSWAP Eq. (5.28), Donelan et al. (1985) Eq. (5.34) and Dobson et al. (1989) Eq. (5.39)."lation results to the theoretical line.Being parameterized by the fetch coordinate, the numerical simulation results evolve from the left to the right on the graph, from the dimensionless fetch equal to 0, to 3.0 × 10 4 .

Comparison with the experiments
A comparison of limited fetch and duration limited simulations with the experimental results by Long and Resio (2007) and the theoretical prediction based on Eq. ( 10) is presented in Figs.11 and 22.One should note that the numerical results and theoretical prediction line with corresponding values of λ fall into the rms deviation (r 2 = 0.939; see Fig. 1) relative to the experimental regression line Eq.( 10).
The dependencies of the dimensionless energy and the frequency on the dimensionless fetch for the limited fetch simulation, superimposed on the experimental observations collected by Young (1999), are presented in Figs.23 and 24, showing good consistency of the presented numerical results and the experimental observations.

Figure 3 .Figure 4 .
Figure 3. Energy local power function index p = d ln E d ln t as a function of dimensionless time tg/U for the wind speed U = 10 m s −1 duration limited case -solid line.Theoretical value of the selfsimilar index p = 10/7 -thick horizontal dashed line.

Figure 5 .Figure 6 .
Figure 5. Mean frequency local power function index −q = d ln ω d ln tas the function of dimensionless time tg/U for the wind speed U = 10 m s −1 duration limited case (solid line).Theoretical value of selfsimilar exponent q = −3/7 -thick horizontal dashed line.

Figure 7 .
Figure 7. Decimal logarithm of the angle averaged spectrum as the function of the decimal logarithm of the frequency for the wind speed U = 10 m s −1 duration limited case -solid line.Spectrum ∼ f −4 -dashed line; spectrum ∼ f −5 -dash-dotted line.

Figure 8 .Figure 9 .
Figure 8. Compensated spectrum as the function of linear frequency f for the wind speed U = 10 m s −1 duration limited case.

Figure 7 Figure 10 .Figure 11 .
Figure 7 presents an angle-integrated energy spectrum as the function of frequency, in logarithmic coordinates.One can see that it consists of the segments of the spectral peak region, the inertial (equilibrium) range ω −4 spanning from the spectral peak to the beginning of the "implicit dissipation" f d = 1.1 Hz, and a Phillips high-frequency tail ω −5 starting approximately from f d = 1.1 Hz.

Figure 14 .
Figure 14.Dimensionless mean frequency as the function of the dimensionless fetch, calculated as f = 1 2π

Figure 15 .Figure 16 .
Figure 15.Local mean frequency exponent −q = d ln ω d ln x as the function of dimensionless fetch xg/U 2 for the limited fetch case.Wind speed U = 10 m s −1 -solid line; wind speed U = 5 m s −1dashed line.Horizontal dashed line -target value of the self-similar exponent q = 0.3.

Figure 17 .Figure 18 .
Figure 17.Decimal logarithm of the angle averaged spectrum as the function of the decimal logarithm of the frequency for the wind speed U = 10 m s −1 limited fetch case -solid line.Spectrum ∼ f −4 -dashed line; spectrum ∼ f −5 -dash-dotted line.

Figure 20 .Figure 21 .Figure 22 .
Figure 20.Compensated spectrum as the function of linear frequency f for the wind speed 10 m s −1 limited fetch case.

Figure 23 .
Figure 23.The solid line (pointed to by arrows) presents non-dimensional total energy from the limited fetch numerical experiments, superimposed onto Fig. 5.4, which is adapted fromYoung (1999).The original caption is "A composite of data from variety of studies showing the development of the non-dimensional energy, ε as a function of non-dimensional fetch, χ .The original JONSWAP study(Hasselmann  et al., 1973)  used the data marked, JONSWAP, together with that ofBurling (1959) and Mitsuyasu (1968).Also shown are a number of growth curves obtained from the various data sets.These include JONSWAP Eq. (5.27),Donelan et al. (1985)  Eq. (5.33) andDobson et al. (1989)  Eq. (5.38)."

Figure 24 .
Figure24.The solid line (pointed to by arrows) presents non-dimensional average frequency as the function of the fetch for limited fetch numerical experiments, superimposed on Fig.5.5 adapted fromYoung (1999).The original caption is "A composite of data from a variety of studies showing the development of the non-dimensional peak frequency, ν as a function of non-dimensional fetch, χ .The original JONSWAP study(Hasselmann et al., 1973)  used all the data shown with the exception of that marked Donelan et al. (1985) and Dobson et al. (1989).Also shown are a number of growth curves obtained from the various data sets.These include JONSWAP Eq. (5.28), Donelan et al. (1985) Eq. (5.34) and Dobson et al. (1989) Eq. (5.39)." Thomson et al. (2013)ced source terms for wave generation however,Thomson et al. (2013)recently used extensive data from Ocean Station Papa to show that there was minimal wind input into the wave spectrum in the equilibrium range.