High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind

Abstract

The paper describes the numerical experiments with a three-dimensional phase-resolving model based on the initial potential equation of motion with free surface at deep water in the periodic domain written in the surface-following nonstationary curvilinear nonorthogonal coordinate system. The numerical scheme is based on Fourier-transform method. The vertical velocity on surface is calculated by solving the three-dimensional Poisson equation for the velocity potential. The velocity potential is represented as a sum of linear and nonlinear components. The linear component is described by Laplace equation. The nonlinear component is calculated by solution of the three-dimensional Poisson equation with the iterated right-hand side. The model includes some algorithms for calculation of the energy input from wind as well as for calculation of breaking and high-frequency dissipation. Initially, the conditions are assigned as a set of small waves corresponding to JONSWAP spectrum at high wave number. In response to waves’ growth, the spectrum shifts to lower wave numbers. The evolution of spectrum is generally in an agreement with the observed data. The wave spectrum and the spectra of different rates of energy transformation as well as the statistical characteristics of wave field for different stages of development are described.

1. Introduction

The development of waves under the action of wind is a process that is difficult to simulate since surface waves are very conservative and their energy changes over hundreds and thousands of periods. This is why the most popular method is spectral modeling based on the averaged over phase equations for spectral energy. In this approach, the waves as physical objects are actually absent since the evolution of spectral distribution of the wave energy is simulated. The description of input and dissipation in this approach is not directly connected with the formulation of the problem, but it is rather adopted from other branches of the wave theory where waves are the objects of investigation. However, the spectral approach was found to be the only method capable to describe the space and time evolution of wave field in the ocean. The phase resolving models (or ‘direct’ models) designed for reproducing waves themselves cannot compete with spectral models since such models typically can reproduce the evolution of just several thousands of large waves. Nevertheless, the direct wave modeling plays an ever-increasing role in the geophysical fluid dynamics because it gives the possibility to investigate the processes that cannot be reproduced by spectral models.

The spectral model assumes that wave field consists of a superposition of linear waves with random phases and arbitrary angle distribution. Being converted to a physical wave field, it looks unreal because real waves usually have prolonged smooth troughs and sharp peaks. Such shape suggests that the waves are similar to Stokes waves. For any given wave spectrum, the wave field can be represented as a superposition of linear Fourier modes with random phases [1]. It can be represented also as a superposition of Stokes modes. The calculations of statistical characteristics for both wave fields show that they are nearly identical. However, such conclusion could be made with no calculations because typical steepness of sea waves at reasonable spectral resolution is of the order of 0.01–0.05, so all the amplitudes of Stokes modes starting from the second one are small. It follows that the specific shape of sea waves is a dynamic property.

The breaking is usually initiated in the vicinity of wave peaks, so the breaking parameterization describes the processes, which are, in principle, impossible in a linear wave field. The breaking is concentrated in the separated narrow intervals; an instantaneous spectrum is discrete and is shifted to high frequencies. The spectral description of the isolated extreme waves is also impossible. The Fourier transform of a large-scale wave field including separated large waves does not provide any indications of their appearance.

The input energy to waves is based on the assumption that each mode induces a pressure mode with a certain amplitude and phase. In reality, the pressure field is quite complicated and not directly connected with the surface elevation because of the systematic separation of air flow behind wave crests (which was shown experimentally in Refs. [2, 3] with the coupled wind-wave model). There are many other complications in the coupled wind-wave dynamics [4] including the processes connected with sprays, bubbles, foam, and structure of the high-frequency wave spectrum.

Most (but not all) of the processes mentioned above can be investigated using the numerical modeling that is a perfect instrument for development of parameterization of physical processes for spectral wave models.

The phase-resolving models can completely replace the spectral models for direct simulation of wave regimes of small water basins, for example, port harbors (see Refs. [5, 6]). Other approaches of direct modeling are discussed in Refs. [7, 8].

Over the past decades, a big volume of papers devoted to the numerical methods developed for investigation of wave processes has been published. The most advanced among them are Finite Difference Method [5, 6], Finite Volume Method [9], Finite Element Method [10, 11], Boundary (Integral) Element Method [12], Spectral HOS Methods [13, 14, 15, 16, 17], the Smoothed Particle Hydrodynamics Method [18], Large Eddy Simulation Method (LES) [19, 20], Moving Particle Semi-Implicit Method [21], Constrained Interpolation Profile Method [22], and Method of Fundamental Solutions [23]. Most of the models were designed for engineering application handling such processes as overturning waves, broken waves, waves generated by landslides, freak waves, solitary waves, tsunamis, violent sloshing waves, interaction of extreme waves with beaches, as well as interaction of steep waves with fixed and different floating structures. The wave models designed for engineering applications seem to be more advanced than the models for pure geophysical research. However, as a rule, the engineering models pay little attention to description of physical processes that are responsible for a long-term evolution of wave spectrum. A more detailed review of direct numerical models is given in Ref. [8].

Until recently, the direct modeling was used for reproduction of a quasi-stationary wave regime when wave spectrum does not change significantly. An example of direct numerical modeling of surface wave evolution is given in Ref. [24] where the development of wave field was calculated by using of a two-dimensional model based on full potential equations written in the conformal coordinates. A model included the algorithms for parameterization of the input and dissipation of energy (a description of similar algorithms is given below). The model successfully reproduced the evolution of wave spectrum under the action of wind. That model was a prototype of 3-D model, because being very fast it was convenient for development of the physical process parameterization. However, the strictly one-dimensional (unidirected) waves are not quite realistic since the unidirected waves in the presence of the small-amplitude perturbations relatively quickly turn into the two-dimensional wave field [25]. Hence, the full problem of wave evolution should be formulated on the basis of three-dimensional equations. Such 3-D calculations were done by Chalikov [26]. The model included parameterization of the main physical processes: input energy, different types of dissipation, and transformation of spectrum due to the nonlinear interaction. The last process does not require any parameterization because the nonlinearity is described with equations. The model used a relatively poor resolution ($1024\times 512$ nodes in x and y directions). However, the calculations reproduced the evolution of wave spectrum and the spectra of the main physical processes such as input and dissipation of energy and nonlinear interactions. As long as we know, it was the first attempt to reproduce the development of waves based on full three-dimensional equations with a direct solution of 3-D equation for the velocity potential. The current paper is devoted to development of the method, including the tuning and modifications of the algorithm, the increase of resolution, and the integration for longer periods. The most important but pure technical modifications were introduced in the numerical scheme for Poisson equation for the velocity potential. The algorithms for calculation of the input and dissipation remain nearly the same, but the numerical parameters in those schemes were changed to achieve a better agreement with the rate of spectrum evolution given by JONSWAP approximation.

2. Equations

The nonstationary surface-following nonorthogonal coordinate system is used as follows:

$$\xi =x,\vartheta =y,\zeta =z-\eta \left(\xi ,\vartheta ,\tau \right),\tau =t,$$E1

where $\eta \left(x,y,t\right)=\eta \left(\xi ,\vartheta ,\tau \right)$ is a moving periodic wave surface given by the Fourier series

$$\eta \left(\xi ,\vartheta ,\tau \right)=\sum _{-M_x<k<M_x}\sum _{-M_y<l<M_y}{h}_{k,l}\left(\tau \right){\mathrm{\Theta }}_{k,l},$$E2

where k and l are the components of the wave number vector $\mathrm{k}$; $h_{k,l}\left(\tau \right)$ are Fourier amplitudes for elevations $\eta \left(\xi ,\vartheta ,\tau \right)$; $M_x$ and $M_y$ are the numbers of modes in the directions $\xi$ and $\vartheta$, respectively, while ${\mathrm{\Theta }}_{k,l}$ are the Fourier expansion basis functions represented as the matrix:

$${\mathrm{\Theta }}_{\mathit{kl}}=\left\{\begin{array}{cc}cos\left(\mathit{k\xi }+\mathit{l\vartheta }\right)& -M_x\le k\le M_x,-M_y<l<0\\ cos\left(\mathit{k\xi }\right)& -M_x\le k\le 0,l=0\\ sin\left(\mathit{k\xi }\right)& 0\le k\le M_y,l=0\\ sin\left(\mathit{k\xi }+\mathit{l\vartheta }\right)& -M_x\le k\le M_x,0<l\le M_y\end{array}\right.$$E3

The 3-D equations of potential waves in the system of coordinates (1) at $\zeta \le 0$ take the following form:

$${\eta }_{\tau }=-{\eta }_{\xi }{\phi }_{\xi }-{\eta }_{\vartheta }{\phi }_{\vartheta }+\left(1+{\eta }_{\xi }^2+{\eta }_{\vartheta }^2\right){\mathrm{\Phi }}_{\varsigma },$$E4

$${\phi }_{\tau }=-\frac{1}{2}\left({\phi }_{\xi }^2+{\phi }_{\vartheta }^2-\left(1+{\eta }_{\xi }^2+{\eta }_{\vartheta }^2\right){\mathrm{\Phi }}_{\zeta }^2\right)-\eta -p,$$E5

$${\mathrm{\Phi }}_{\mathit{\xi \xi }}+{\mathrm{\Phi }}_{\mathit{\vartheta \vartheta }}+{\mathrm{\Phi }}_{\mathit{\zeta \zeta }}=\mathrm{\Upsilon }\left(\mathrm{\Phi }\right),$$E6

where $\mathrm{\Upsilon }$ is the operator:

$$\mathrm{\Upsilon }\left(\right)=2{\eta }_{\mathit{\xi }}{\left(\right)}_{\mathit{\xi \zeta }}+2{\eta }_{\mathit{\vartheta }}{\left(\right)}_{\mathit{\vartheta \zeta }}+\left({\eta }_{\mathit{\xi \xi }}+{\eta }_{\mathit{\vartheta \vartheta }}\right){\left(\right)}_{\zeta }-\left({\eta }_{\xi }^2+{\eta }_{\vartheta }^2\right){\left(\right)}_{\mathit{\zeta \zeta }},$$E7

capital fonts $\mathrm{\Phi }$ are used for the domain $\zeta <0$, while the lower case $\phi$ refers to $\zeta =0.$ The term p in Eq. (5) describes the pressure on the surface $\zeta =0$.

It is suggested by Chalikov et al. [7] that it is convenient to represent the velocity potential $\phi$ as a sum of the two components: a linear component $\overline{\mathrm{\Phi }},\left(\overline{\phi }=\overline{\mathrm{\Phi }}\left(\xi ,\vartheta ,0\right)\right)$ and an arbitrary nonlinear component $\tilde{\mathrm{\Phi }},\left(\tilde{\phi }=\tilde{\mathrm{\Phi }}\left(\xi ,\vartheta ,0\right)\right)$:

$$\phi =\overline{\phi }+\tilde{\phi },\mathrm{\Phi }=\overline{\mathrm{\Phi }}+\tilde{\mathrm{\Phi }}.$$E8

The linear component $\overline{\mathrm{\Phi }}$ satisfies Laplace equation:

$${\overline{\mathrm{\Phi }}}_{\mathit{\xi \xi }}+{\overline{\mathrm{\Phi }}}_{\mathit{\vartheta \vartheta }}+{\overline{\mathrm{\Phi }}}_{\mathit{\zeta \zeta }}=0,$$E9

with a known solution:

$$\overline{\mathrm{\Phi }}\left(\xi ,\vartheta ,\zeta ,\tau \right)=\sum _{k,l}{\overline{\phi }}_{k,l}\left(\tau \right)exp\left(\left|k\right|\zeta \right){\mathrm{\Theta }}_{k,l},$$E10

where $\left|k\right|={\left(k^2+l^2\right)}^{1/2}$, ${\overline{\phi }}_{k,l}$ are the Fourier coefficients of the surface linear potential $\overline{\phi }$ at $\zeta =0).$ The solution satisfies the following boundary conditions:

$$\begin{array}{l}\varsigma =0:\overline{\mathrm{\Phi }}=\overline{\phi }\\ \varsigma \to -\infty :{\tilde{\mathrm{\Phi }}}_{\zeta }\to 0\end{array}$$E11

The nonlinear component satisfies the equation:

$${\tilde{\mathrm{\Phi }}}_{\mathit{\xi \xi }}+{\tilde{\mathrm{\Phi }}}_{\mathit{\vartheta \vartheta }}+{\tilde{\mathrm{\Phi }}}_{\mathit{\zeta \zeta }}=\mathrm{\Upsilon }\left(\tilde{\mathrm{\Phi }}\right)+\mathrm{\Upsilon }\left(\overline{\mathrm{\Phi }}\right),$$E12

Eq. (12) is solved with the boundary conditions:

$$\begin{array}{l}\varsigma =0:\tilde{\mathrm{\Phi }}=0\\ \varsigma \to -\infty :{\tilde{\mathrm{\Phi }}}_{\zeta }\to 0\end{array}$$E13

The presentation (8) is not used for solution of the evolutionary Eqs. (4) and (5) because it does not provide any improvements of accuracy and speed.

Eqs. (4)–(6) are written in a nondimensional form by using the following scales: length L, where 2πL is a (dimensional) period in the horizontal direction; time L^1/2g^−1/2; and velocity potential $L^{3/2}{g}^{1/2}$ (g is the acceleration of gravity). The pressure is normalized by the water density, so that the pressure scale is Lg. Eqs. (4)–(6) are self-similar to the transformation with respect to L. The dimensional size of the domain is $2\mathit{\pi L}$, so the scaled size is $2\pi .$ All of the results presented in this paper are nondimensional. Note that the number of the Fourier modes can be different in the x and y directions. In this case, it is assumed that the two-length scales $L_x$ and $L_y$ are used. The nondimensional length of the domain in the y-direction remains equal to $2\pi$, and a factor $r=L_x/L_y$ is introduced into the definition of a differential operator in the Fourier space.

The derivatives of a linear component $\overline{\mathrm{\Phi }}$ in (7) are calculated analytically. The scheme combines a 2-D Fourier transform method in the ‘horizontal surfaces’ and a second-order finite-difference approximation on the stretched staggered grid defined by the relation $\mathrm{\Delta }{\zeta }_{j+1}=\chi \mathrm{\Delta }{\zeta }_j$ ($\mathrm{\Delta }\zeta$ is a vertical step, while $j=1$ at the surface). The stretched grid provides an increase in accuracy of approximation for the exponentially decaying modes. The values of the stretching coefficient $\chi$ lie for different settings within the interval 1.01–1.20; in the current work, the value$\gamma =1.2$ was used at the number of levels $L_w=10$. Such poor resolution was possible to use because of the separation of the potential into a large linear and a small nonlinear part, so Eq. (12) was used only for calculation of a small correction for the potential. A high value of the stretching coefficient provided high resolution in the vicinity of surface for accurate calculations of the surface vertical derivative for the potential. The finite-difference second-order approximation of vertical operators in Eq. (12) on a nonuniform vertical grid is quite straightforward (see Ref. [8]). Eq. (12) is solved as Poisson equations with the iterations over the right-hand side by TDMA method [27]. At each time step, the iterations start with a right-hand side calculated at the previous time step. A relative accuracy of the solution in terms of the vertical derivative of the potential on the surface was equal to ${10}^{-6}$. The typical number of iterations was 2–5.

The accuracy of the adiabatic version of equations was validated by reproducing a moving Stokes wave with the steepness $\mathit{AK}=0.40$ ($А$ is a half of trough-to-crest wave height; $K=1$ is the wave number of the first mode). An algorithm for calculation of Stokes wave with the prescribed accuracy was suggested by Chalikov and Sheinin [28]. The scheme based on the conformal coordinates is very effective: the calculations were carried out in 100 ms at notebook (2.10 GHz). The dependence of Stokes wave spectra on the wave number is shown in Figure 1. In fact, about 2000 curves obtained in the course of calculations, with the interval $\mathrm{\Delta }t=1$, were plotted. Due to improvement of the numerical scheme, the accuracy of reproduction of Stokes wave is considerably higher than for the scheme used in Refs. [7, 8].

Figure 1.

The 2000 spectra of Stokes wave with steepness $\mathit{AK}=0.40$ as a function of wave number $k$ superimposed on each other. The gray area in the right left corner is a computational noise.

As seen, up to $S\sim {10}^{-12}\left(k=22\right)$, the spectra of Stokes waves remain with high accuracy the same as it was assigned in initial conditions. At higher frequencies, the random disturbances appear. Note that this validation is not trivial: even small inaccuracies in the numerical scheme cause a fast distortion of the spectra, like in the bottom part of Figure 1. We consider these results as a serious evidence of high accuracy of the adiabatic version of the model. The previous version of 3-D model [26] allowed carrying out a long simulation of Stokes wave not steeper than $\mathit{AK}=0.30$. The right-hand sides of Eqs. (4) and (5) were calculated with a use of Fourier transform method: the nonlinear terms were calculated at the extended grid with size $\left(4M_x\times 4M_y\right)$, and then by the inverse Fourier transform, they were returned to the Fourier grid. The fourth-order Runge-Kutta scheme was used for integration in time. The equation for the potential was solved at each of the four substeps of time step.

The simulations described by Chalikov [26] were a first attempt to reproduce the development of wave field assigned in the initial conditions as a group of small waves at high wave number under the action of strong wind. The initial elevation was generated as a superposition of linear waves corresponding to JONSWAP spectrum [29] with random phases. The initial Fourier amplitudes for the surface potential were calculated by the formulas of the linear wave theory. The details of the initial conditions are of no importance because the initial energy level is quite low. The wave peak was placed to the wave number equal to 100. The wind velocity was assigned equal to$4c_{100}$, where $c_{100}$ is a phase velocity of the 100th mode. A detailed description of the scheme and its validation is given in Refs. [7, 8].

The simulation described in the current paper was performed with a doubled resolution in both directions, with the improved numerical scheme for Poisson equations and modified parameters in the scheme for calculations of energy transitions.

3. Energy input

The detailed description of the algorithm for calculation of energy input is given in Ref. [26]. The energy and momentum are transferred from air to water by the surface pressure field and tangent stress. According to the most reasonable theory [30], the Fourier components of surface pressure $p$ are connected with those of the surface elevation through the following expression:

$$p_{k,l}+\mathrm{i}{p}_{-k,-l}=\frac{{\rho }_a}{{\rho }_w}\left({\beta }_{k,l}+\mathrm{i}{\beta }_{-k,-l}\right)\left(h_{k,l}+\mathrm{i}{h}_{-k,-l}\right),$$E14

where $h_{k,l},h_{-k,-l},{\beta }_{k,l},{\beta }_{-k,-l}$ are real and imaginary parts of elevation $\eta$, and the so-called $\beta$-function, ${\rho }_a/{\rho }_w$, is the ratio of air and water densities. Both $\beta$ coefficients are the functions of the nondimensional frequency

$$\mathrm{\Omega }=\omega U/g,$$E15

that characterizes the ratio of wind velocity to phase velocity of $c_k$:

$$\mathrm{\Omega }=U/c_k$$E16

Since the supplying of wave with the energy and momentum occurs in a layer whose height is proportional to the wave length, it is reasonable to suggest that the reference height for the wind velocity should be different for a different virtual wave length (distance ${\lambda }_k/cos{\theta }_i$ between the wave peaks in wind direction; the index $i$ denotes a direction of mode). The wind velocity can be found by interpolation or extrapolation to the level:

$$z_{i,k}=0.5{\lambda }_k/cos{\theta }_i$$E17

The definition of ${\mathrm{\Omega }}_{к}$ should take into account the angle ${\theta }_i$ between the vector $U$ and the direction of wave mode. Finally, the virtual nondimensional frequency takes the form:

$${\mathrm{\Omega }}_{i,k}={\omega }_{k}cos{\theta }_{i}U\left(z_k\right)/g=cos{\theta }_{i}U\left(z_k\right)/c_k$$E18

where $c_k=g/{\omega }_k$ is the phase velocity of kth mode.

For experimental derivation of the shape of β-function, it is necessary to simultaneously measure the wave surface elevation and nonstatic pressure on the surface [31, 32, 33, 34, 35]. The data obtained in this way allowed constructing an imaginary part of β-function used in some versions of the wave forecasting models [36]. The data on experimental $\beta$-function are compared in Ref. [4]. The values of $\beta$ within the interval $\left(0<\mathrm{\Omega }<10\right)$ differ by decimal orders. Hence, the question arises: in what way, using such a different input, the spectral models provide a reasonable agreement with the observations. The answer is very simple: the researchers have the possibility to modify the parameterization of dissipation. Despite the hundreds of papers, the knowledge on dissipation is even poorer than the knowledge on the energy input. Finally, only the sum of those source terms regulates the growth of total wave energy. Such situation is far from being perfect since the energy input and dissipation have totally different spectral properties.

The second way of the β-function evaluation is based on the results of numerical investigations of the statistical structure of the boundary layer above waves with the use of Reynolds equations and an appropriate closure scheme. In general, this method works so well that many problems in the technical fluid mechanics are often solved not experimentally but by using the numerical models [37, 38]. This method was being developed beginning from Refs. [39, 40] and followed by Refs. [41, 42, 43]. The results were implemented in the WAVEWATCH model, i.e., the third-generation wave forecast model [44], and thoroughly validated against the experimental data in the course of developing WAVEWATCH-III [45]. Most of the schemes for the calculations of $\beta$-function consider a relatively narrow interval of the nondimensional frequencies $\mathrm{\Omega }$. In the current work, the range of frequencies covers the interval $\left(0<\mathrm{\Omega }<10\right)$, and occasionally, the values of $\mathrm{\Omega }>10$ can appear.

The most reliable data on $\beta$-function are concentrated in the interval $-10<\mathrm{\Omega }<10$ (the negative values of $\mathrm{\Omega }$ correspond to the wave modes running against wind). In the current calculations, the modes running against wind are absent. The function $\beta$ can be approximated by the formulas:

$${\beta }_i=\left\{\begin{array}{cc}{\beta }_0+a_0\left(\mathrm{\Omega }-{\mathrm{\Omega }}_0\right)+a_1{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_0\right)}^2& \mathrm{\Omega }>{\mathrm{\Omega }}_0\\ {\beta }_0-a_0\left(\mathrm{\Omega }-{\mathrm{\Omega }}_0\right)+a_1{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_0\right)}^2& \mathrm{\Omega }<{\mathrm{\Omega }}_0\end{array}\right.,$$E19

$${\beta }_r=\left\{\begin{array}{cc}{\beta }_1+a_3\left(\mathrm{\Omega }-{\mathrm{\Omega }}_2\right)& \mathrm{\Omega }<{\mathrm{\Omega }}_2\\ a_2{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_1\right)}^2& {\mathrm{\Omega }}_2<\mathrm{\Omega }<{\mathrm{\Omega }}_3\\ {\beta }_1-a_3\left(\mathrm{\Omega }-{\mathrm{\Omega }}_3\right)& \mathrm{\Omega }>{\mathrm{\Omega }}_3\end{array}\right.,$$E20

where ${\mathrm{\Omega }}_0=0.355$, ${\mathrm{\Omega }}_1=1.20$, ${\mathrm{\Omega }}_2=-18.8$, ${\mathrm{\Omega }}_3=21.2$, $a_0=0.0228,$$a_1=0.0948,$$a_2=-0.372,$$a_3=14.8,$${\beta }_0=-0.02,$${\beta }_1=-148.0$.

The wind velocity remains constant throughout the integration. The values of $\mathrm{\Omega }$ for other wave numbers are calculated by assuming that the wind profile is logarithmic.

Note that the formulation of wind and waves interaction can be significantly improved by coupling the wave model with the 1-D Wave Boundary Layer model [4]. The next step can be the coupling of wave model with the 3-D model of WBL based on the closure schemes or LES model (see Ref. [46]).

4. Energy dissipation

The current version of the model includes three types of dissipation (see details in Ref. [26]).

1. The energy can decrease due to the errors of approximations in space and time that depend on the number of Fourier modes, number of knots in the physical space, the vertical grid used for approximation of Poisson equation (6), and the criterion for accuracy of its solution. All of those errors that produce the ‘numerical dissipation’ can be referred to the adiabatic part of the models (4)–(6) at $p=0$. The rate of this dissipation can be reduced by the use of a better resolution and a higher accuracy of approximation, but this way leads to deceleration of the calculations with the model already running for a very long time.

   Opposite to the numerical dissipation, there exists another type of energy loss that has rather a physical nature. The nonlinear interaction of different modes forms a flux of energy directed outside of the computational domain. We call it the ‘nonlinear dissipation.’ The numerical and nonlinear dissipation can hardly be considered separately. The estimation of rate of the numerical/nonlinear dissipation can be easily done by the comparison of full energy before and after the time step for the adiabatic part of the model (see Section 4 in Ref. [26]). In the current calculations, the loss of energy for one time step was about ${10}^{-4}\%$, which is by 2–3 orders less that the rate of energy change due to input energy. Since we prefer to consider the process described by Eqs. (4)–(6) as adiabatic one, at each time step we restore the energy lost by both the numerical and nonlinear dissipation.

2. A long-term integration of full fluid mechanics equations always shows the spreading of spectrum to both high and low frequencies (wave numbers). The nonlinear flux of energy directed to the small wave numbers produces downshifting of spectrum, while an opposite flux forms a shape of the spectral tail. The second process that we call the ‘tail dissipation’ can produce accumulation of energy near the ‘cut’ wave number. The growth of amplitudes at high wave numbers is followed by growth of the local steepness and development of the numerical instability. To support the stability, additional terms are included into the right-hand sides of Eqs. (4) and (5):

   $$\frac{\partial {\eta }_{k,l}}{\partial \tau }=E_{k,l}-{\mu }_{k,l}{\eta }_{k,l},$$E21
   $$\frac{\partial {\phi }_{k,\mathrm{l}}}{\partial \tau }=F_{k,l}-{\mu }_{k,l}{\phi }_{k,l}$$E22

   (where $E_{k,l}$ and $F_{k,l}$ are the Fourier amplitudes of the right-hand sides of Eqs. (4) and (5); the value of ${\mu }_{k,l}$ is equal to zero inside the ellipse with semi-axes $d_m{M}_x$ and $d_m{М}_y$; then, it grows quadratically with $\left|k\right|$ up to the value $c_m$ and is equal to $c_m$ outside of the outer ellipse (see details in Ref. [26]). This method of filtration that we call the ‘tail dissipation’ was developed and validated with the conformal model [28]. The sensitivity of the results to the parameters in Refs. (21) and (23) is not large. The aim of the algorithm is to support smoothness and monotonicity of the wave spectrum within the high wave number range.

3. The main process of wave dissipation is the ‘breaking dissipation.’ This process is taken into account in all the spectral wave forecasting models similar to WAVEWATCH (see Refs. [44, 47]). Since there are no waves in the spectral models, no local criteria of wave breaking can be formulated. This is why the breaking dissipation is represented in the spectral models in a distorted form. The real breaking occurs in the relatively narrow areas of the physical space; however, the spectral image of such breaking is stretched over the entire wave spectrum, while in reality, the breaking decreases height and energy of separate waves. This contradiction occurs because the waves in the spectral models are assumed to be linear. In fact, a nonlinear sharp wave breaks in the physical space. Such wave is often composed of several local modes. It is clear that the state-of-art wave models should account for the threshold behavior of a breaking wave, that is, waves will not break unless their steepness exceeds the threshold [48, 49, 50].

The instability of the interface leading to breaking is an important though poorly developed problem of fluid mechanics. In general, this essentially nonlinear process should be investigated for the two-phase flow. Such approach was demonstrated, for example, by Iafrati [51].

The problem of breaking parameterization includes two points: (1) establishment of a criterion of the breaking onset and (2) development of the algorithm of the breaking parameterization. The problem of breaking is discussed in details in Ref. [47]. It was found in Ref. [52] that the clear predictor of breaking formulated in dynamical and geometrical terms, probably, does not exist. The consideration of the exact criterion for the breaking onset for the models using transformation of the coordinate type of (1) is useless since the numerical instability in such models occurs not because of the approach of breaking but because of the appearance of the high local steepness. The description of breaking in the direct wave modeling should satisfy the following conditions: (1) it should prevent the onset of instability at each point of millions of grid points over many thousands of time steps; (2) it should describe in a more or less realistic way the loss of the kinetic and potential energies with preservation of balance between them; and (3) it should preserve the volume. It was suggested by Chalikov [53] that an acceptable scheme can be based on a local highly selective diffusion operator with a special diffusion coefficient. Several schemes of such type were validated, and finally, the following scheme was chosen:

$${\eta }_{\tau }=E_{\eta }+J^{-1}\left(\frac{\partial }{\partial \xi }{B}_{\xi }\frac{\partial \eta }{\partial \xi }+\frac{\partial }{\partial \vartheta }{B}_{\vartheta }\frac{\partial \eta }{\partial \vartheta }\right),$$E23

$${\phi }_{\tau }=F_{\phi }+J^{-1}\left(\frac{\partial }{\partial \xi }{B}_{\xi }\frac{\partial \phi }{\partial \xi }+\frac{\partial }{\partial \vartheta }{B}_{\vartheta }\frac{\partial \phi }{\partial \vartheta }\right),$$E24

where $F_{\eta }$ and $F_{\phi }$ are the right-hand sides of Eqs. (4) and (5) including the tail dissipation terms; $B_{\xi }$ and $B_{\vartheta }$ are the diffusion coefficients. The probability of high negative values of the curvilinearity is by orders larger than the probability calculated over the ensemble of linear modes with the spectra generated by the nonlinear model.

The curvilinearity turned out to be very sensitive to the shape of surface. This is why it was chosen as a criterion of the approaching breaking. The coefficients $B_{\xi }$ and $B_{\vartheta }$ depend nonlinearly on the curvilinearity

$$B_{\xi }=\left\{\begin{array}{cc}{C}_B{\eta }_{\mathit{\xi \xi }}^2& {\eta }_{\mathit{\xi \xi }}<{\eta }_{\mathit{\xi \xi }}^{\mathit{cr}}\\ 0& {\eta }_{\mathit{\xi \xi }}\ge {\eta }_{\mathit{\xi \xi }}^{\mathit{cr}}\end{array}\right.$$E25

$$B_{\vartheta }=\left\{\begin{array}{cc}{C}_B{\eta }_{\mathit{\vartheta \vartheta }}^2& {\eta }_{\mathit{\vartheta \vartheta }}<{\eta }_{\mathit{\xi \xi }}^{\mathit{cr}}\\ 0& {\eta }_{\mathit{\vartheta \vartheta }}\ge {\eta }_{\mathit{\xi \xi }}^{\mathit{cr}}\end{array}\right.$$E26

where the coefficients at $C_B=0.05$, ${\eta }_{\mathit{\xi \xi }}^{\mathit{cr}}={\eta }_{\mathit{\vartheta \vartheta }}^{\mathit{cr}}=-50$. The algorithm (24)–(27) does not change the volume and decreases the local potential and kinetic energies. It is assumed that the lost momentum and energy are transferred to the current and turbulence (see Ref. [42]). Besides, the energy also goes to other wave modes. The choice of parameters in Refs. (24)–(27) is based on simple considerations: the local piece of surface can closely approach the critical curvilinearity but not exceed it. The values of the coefficients were chosen in the course of multiple experiments to provide agreement with the rate of spectrum development given by JONSWAP approximation.

5. Evolution of wave field

The integration was done for 1,200,000 steps with the time step $\mathrm{\Delta }t=0.005$ up to the nondimensional time $T=6000$, which corresponded to 9550 initial wave peak periods. The total energy of wave motion $Е=E_p+E_k$ ($E_p$ is the potential energy, while $E_k$ is the kinetic energy) is calculated with the following formulas:

$$E_p=0.25\overline{{\eta }^2},E_k=0.5\overline{\overline{\left({\phi }_x^2+{\phi }_y^2+{\phi }_z^2\right)}},E=E_p+E_k,$$E27

where a single bar denotes the averaging over the $\xi$ and $\vartheta$ coordinates, while a double bar denotes the averaging over the entire volume. The derivatives in Ref. (27) are calculated according to the transformation (1). An equation of the integral energy E evolution can be represented in the following form:

$$\frac{\mathit{dЕ}}{\mathit{dt}}=\overline{I}+\overline{D_b}+\overline{D_t}+\overline{N},$$E28

where $\overline{I}$ is the integral input of energy from wind (Eqs. (14)–(20)); $\overline{D_b}$ is a rate of the energy dissipation due to wave breaking (Eqs. (23)–(26)); $\overline{D_t}$ is a rate of the energy dissipation due to filtration of high-wave number modes (‘tail dissipation,’ Eqs. (21) and (22)); $\overline{N}$ is the integral effect of the nonlinear interactions described by the right-hand side of the equations when the surface pressure $p$ is equal to zero. The differential forms for calculation of the energy transformations can be, in principle, derived from Eqs. (4)–(6), but here a more convenient and simple method was applied. Different rates of the integral energy transformations can be calculated with the help of fictitious time steps (i.e., apart from the basic calculations). For example, the value of $\overline{I}$ is calculated by the following relation:

$$\overline{\overline{I}}=\frac{1}{\mathrm{\Delta }t}\left(\overline{\overline{E^{t+\mathrm{\Delta }t}}}-\overline{\overline{E^t}}\right),$$E29

where $\overline{\overline{E^{t+\mathrm{\Delta }t}}}$ is the integral energy of a wave field obtained after one time step with the right side of Eq. (6) containing only the surface pressure calculated with Eqs. (14)–(18).

The evolution of the characteristics calculated by formula (29) is shown in Figure 2.

Figure 2.

The evolution of integral characteristics. The rate of evolution of the integral energy multiplied by ${10}^8$ due to: (1) nonlinear interaction $\overline{\overline{I}}$ (Eq. (29)); (2) tail dissipation ${\overline{\overline{D}}}_t$(Eqs. (21) and (22)); (3) breaking dissipation $D_b$ (Eqs. (23)–(26)); (4) input of energy from wind $\overline{\overline{I}}$ (Eqs. (14)–(20)); and (5) balance of energy $\overline{\overline{I}}+{\overline{\overline{D}}}_t+{\overline{\overline{D}}}_b$. Curve 6 shows the evolution of wave energy ${10}^{5}E$. Vertical gray bars show instantaneous values; thick curve shows the smoothed behavior.

The sharp variation of all the characteristics at $t<500$ is explained by adjustment of the linear initial fields to the nonlinearity. The integral effect of the nonlinear interaction $\overline{\overline{I}}$ (straight line 1) was very close to zero. The tail dissipation $\overline{\overline{D_t}}$ (curve 2) is smaller than the breaking dissipation $\overline{\overline{D_t}}$(curve 3). The value of $\overline{\overline{D_b}}$ has significant fluctuations due to introduction of the criteria (25) and (26). The dissipation ${\overline{\overline{D}}}_b+{\overline{\overline{D}}}_t$ absorbs nearly all of the incoming energy, and just a small part of it is going for growth of waves. The balance of energy $\overline{\overline{B}}=\overline{\overline{I}}+{\overline{\overline{D}}}_t+{\overline{\overline{D}}}_b$ (curve 5) fluctuates and approaches zero when energy $\overline{\overline{E}}$ (curve 6 in Figure 2) approaches saturation.

The time evolution of the integral spectral characteristics is presented in Figure 3.

Figure 3.

Dependence of integral characteristics on fetch (Eq. (32)): (1) weighted mean frequency ${\omega }_w$ (Eq. (31)); (2) peak frequency ${\omega }_p$; (3) energy E (Eq. (27)); and (4) approximation of the peak frequency evolution (Eq. (32)).

Curve 1 corresponds to the weighted frequency ${\omega }_w$

$${\omega }_w={\left(\frac{\int \mathit{\text{ωSdkdl}}}{\int \mathit{\text{Sdkdl}}}\right)}^{1/2},$$E30

where integrals are taken over the entire Fourier domain. The value ${\omega }_w$ is not sensitive to the details of spectrum; hence, it well characterizes the position of spectrum and spectral peak shifting. Curve 2 describes the evolution of the spectral maximum. The step shape of curve corresponds to the fundamental property of downshifting. Opposite to common views, the development of spectrum occurs not monotonically, but by appearance of a new maximum at a lower wave number as well as by attenuation of the previous maximum. It is interesting to note that the same phenomenon is also observed in the spectral model [36].

The value of fetch in the periodic problem can be calculated by integration of the peak phase velocity $c_p={\left|k\right|}^{-1/2}$ over time.

$$F={\int }_{t_0}^t{c}_p\mathit{dt}.$$E31

The numerical experiment reproduces the case when development of wave field occurs under the action of a permanent and uniform wind. This case corresponds to the JONSWAP experiment [29]. It is suggested that the frequency of spectral peak changes as $F^{-1/3}$, while the full energy grows linearly with $F.$ Neither of the dependences can be exact since they do not take into account approaching a stationary regime. Besides, the dependence of frequency on fetch is singular at $F=0.$ A more accurate is the approximation:

$${\omega }_p=\frac{75.6}{5.63+F^{1/3}}.$$E32

Obviously, the dependence ${\omega }_p\sim F^{-1/3}$ is valid in a narrower interval of F. As seen, contrary to ${\omega }_w$, the peak frequency changes not monotonically, but by appearance of a new maximum at a lower wave number as well as by attenuation of the previous maximum. It is interesting to note that the same phenomenon is also observed in the spectral model (16). The dependence of the total energy E on fetch F does not look like a linear one, but it is worth to note that the JONSWAP dependence is evidently inapplicable to a very small and large fetch.

On the whole, the evolutions of integral characteristics of the solution shown in Figures 2 and 3 are smoother than those calculated by Chalikov [26]. It can be explained by multiple technical improvements of the numerical scheme and higher resolution.

The evolution of wave spectrum is shown in Figure 4.

Figure 4.

The wave spectra $S_h\left(r\right)$ integrated over angle $\psi$ in the polar coordinates and averaged over the consequent intervals of length about 500 units of the nondimensional time t (thin curves). The spectra are growing and shifting from right to left. Thick dashed curve is the dependence $S\sim {\omega }^{-4}$; dotted curve corresponds to $S\sim {\omega }^{-5}$ .

The 2-D wave spectrum $S\left(k,l\right)$$\left(0\le k\le M_x,-M_y\le l\le M_y\right)$ averaged over nine time intervals of length equal to $\mathrm{\Delta }t\approx 500$ was transferred to the polar coordinates $S_p\left(\psi ,r\right)\left(-\pi /2\le \psi \le \pi /2,0\le r\le M_x\right)$ and then averaged over the angle $\psi$ to obtain the 1-D spectrum$S_h\left(r\right)$:

$$S_h\left(r\right)=\sum S_p\left(\psi ,r\right)r\mathrm{\Delta }\psi .$$E33

The angle $\psi =0$ coincides with the direction of wind $U$, $\mathrm{\Delta }\psi =\pi /180.$ Even the averaged over angle spectrum looks quite irregular and contains multiple holes and peaks. The spectra are smoothed.

The two-dimensional wave spectra are shown in Figure 5, where the ${log}_{10}\left(S\left(\psi ,r\right)\right)$ averaged over the successive eight periods of length $\mathrm{\Delta }t=500$ is given.

Figure 5.

Sequence of 2-D images of $\mathrm{lо}{\mathrm{g}}_{10}\left(S\left(r,\psi \right)\right)$ averaged over the consequent eight periods of length $\mathrm{\Delta }t=500$. The numbers indicate the period of averaging (the first panel marked 0 refers to the initial conditions). The spectra are normalized by the maximum value of spectrum 8. The horizontal axis corresponds to the wave numbers $r=\left|k\right|$; the angles are shown by rays.

The first panel with a mark 0 refers to the initial conditions. The pictures well characterize the downshifting and angle spreading of spectrum due to the nonlinear interactions.

As seen, each spectrum consists of separated peaks and holes. This phenomenon was observed and discussed by Chalikov et al. [7]. The same results were obtained by Chalikov [26]. The repeated calculations with different resolutions showed that such structure of 2-D spectrum is typical. The locations of peaks cannot be explained by the fixed combination of interacting modes, since in different runs (with the same initial conditions but a different set of phases for the modes), the peaks are located in different locations in the Fourier space.

It is interesting to note that while increasing resolution, the patches with low energy extend. It can be supposed that the current and higher resolutions are excessive, and the process can be simulated with a lower resolution. This statement may be too optimistic, but it can be supported by the following arguments. The multi-mode wave mechanics is different from the multi-scale turbulent motion. The modeling of turbulence at increase of resolution just allows reproducing more details of motion. The increase of resolution in a wave model introduces other wave modes with different phase velocities. Due to dispersion, the solutions (i.e., the evolution of surface) in these two cases will be completely different. It means that the solution does not converge with increase of resolution, which makes no sense.

The situation can be saved if upon reaching the optimal resolution, the new added positions for the modes will not obtain the energy and not participate in solution. The existence of such effect should be carefully validated with the exact wave model. If this effect does not exist, it means that the results of simulations depend completely on the resolution, the reliable simulation of individual evolution of wave field being, in fact, impossible.

The method of calculation of the simulated one-dimensional input and dissipation spectra was described by Chalikov [26]; still, it will be explained here once again, though briefly.

The evolution of the integrated over angle $\psi$ wave spectrum $S_h\left(r\right)$ can be described with the equation:

$$\frac{{\mathit{dS}}_h\left(r\right)}{\mathit{dt}}=I\left(r\right)+D_t\left(r\right)+D_b\left(r\right)+N\left(r\right),$$E34

where $I\left(r\right),D_t\left(r\right),D_b\left(r\right)\text{and}N\left(r\right)$ are the spectra of the input energy, tail dissipation, breaking dissipation, and the rate of nonlinear interactions. All of the spectra shown below were obtained by transformation of the 2-D spectra into a polar coordinate $\left(\psi ,r\right)$ and then integrated over the angles $\psi$ within the interval $\left(-\pi /2,\pi /2\right)$. The spectra can be calculated using an algorithm similar to Eq. (29) for integral characteristics. For example, the spectrum of the energy input $I\left(k,l\right)$ is calculated as follows:

$$I\left(k,l\right)=\left(S_c^{t+\mathrm{\Delta }t}\left(k,l\right)-S_c^t\left(k,l\right)\right)/\mathrm{\Delta }t,$$E35

where $S_c\left(k_x,k_y\right)$ is a spectrum of the columnar energy calculated by the relation:

$$S_c\left(k,l\right)=\frac{1}{2}\left(h_{k,l}^2+h_{-k,-l}^2+\underset{-H}{\overset{0}{\int }}\left(u_{k,l}^2+u_{-k,-l}^2+v_{k,l}^2+v_{-k,-l}^2+w_{k,l}^2+w_{-k,-l}^2\right)\mathit{d\zeta }\right)$$E36

where the grid values of velocity components $u,v,w$ are calculated by the relations:

$$u={\phi }_{\xi }+{\phi }_{\zeta }{\eta }_{\xi },v={\phi }_{\vartheta }+{\phi }_{\zeta }{\eta }_{\vartheta },w={\phi }_{\zeta },$$E37

and $u_{k,l},v_{k,l}\text{and}{w}_{k,l}$ are the real Fourier coefficients, while for the negative indices—the imaginary ones.

For calculation of $I\left(k,l\right)$, the fictitious time steps $\mathrm{\Delta }t$ are made only with a term responsible for the energy input, that is, the surface pressure $p$. The spectrum $I\left(k,l\right)$ was averaged over the periods $\mathrm{\Delta }t\approx 500$, then transformed into a polar coordinate system and integrated in the Fourier space over the angles $\psi$ within the interval $\left(-\pi /2,\pi /2\right).$ Such procedure was used for calculation of all the terms in the right side of Eq. (34). In the current version of the model, the calculations of integral (28) and spectral (34) transformations were combined with the calculations of the right sides of Eqs. (4) and (5).

The rates of transformation of spectrum are shown in Figure 6. The integral term describing the nonlinear interaction $\overline{N}$ in Eq. (28) is small (as compared with the local values of $N_{k,l}$), but the magnitude of spectrum $N\left(r\right)$ is comparable with the input $I\left(r\right)$ and dissipation $D_t\left(r\right)$ and $D_b\left(r\right)$ terms (panel 1 in Figure 5). The shape of spectrum $N\left(r\right)$ confirms prediction of the quasi-linear theory [54, 55]. At the low wave number slope of the spectrum, the nonlinear influx of energy is positive, while at the opposite slope, it is negative. This process produces shifting of spectrum to the lower wave number (downshifting). The input of energy due to the nonlinear interactions is observed in a high frequency part of spectrum, which also agrees with Hasselmann’s theory. Note that the nonlinear interactions also produce widening of spectrum.

Figure 6.

The rates of transformation spectra multiplied by ${10}^9$: (1) nonlinear interaction $N\left(r\right)$(1); (2) input energy $I\left(r\right)$; (3) breaking dissipation $D_b\left(r\right)$; (4) tail dissipations $D_t\left(r\right)$ multiplied by ${10}^{11}$. All spectra are obtained by transformation of the 2-D spectra into the polar coordinate $\left(\psi ,r\right)$ and then integrated over the angles $\psi$ within the interval $\left(-\pi /2,\pi /2\right)$.

The spectral distribution of the energy input from wind $I\left(r\right)$ (panel 2 in Figure 6) is in general similar to wave spectrum since it depends linearly on the spectral density (Figure 3). The dissipation rate $D_b\left(r\right)$ is negative (panel 3), and its minimum is shifted a little to higher frequencies from the wave spectrum peak. The tail dissipation (Panel 4) is smaller by two orders than the other terms, but it plays an important role of supporting numerical stability.

The residual rate of transformation of spectrum ${\mathit{dS}}_h\left(r\right)/\mathit{dt}$ averaged over eight consequent periods is shown in Figure 7. The numbers in the top part of panel indicate the averaged wave number of the spectral peak. The second set of numbers refers to the corresponding spectrum of the residual input of energy. As seen, the maximum of the input energy is located to the left of the spectral peak, that is, on a low-wavenumber spectral slope. The obtained energy causes downshifting of spectrum and supports the shape of a high wavenumber slope and spectral tail. In the equilibrium regime, all the incoming energy are consumed for supporting the shape of the entire spectrum.

Figure 7.

The spectral distribution of the residual input of energy ${\mathit{dS}}_h\left(r\right)/\mathit{dt}$ in the energy-containing part of spectrum for different stages of wave development. The numbered vertical lines indicate positions of the spectral peak. The second set of numbers shows the peaks of the corresponding residual input.

However, the dynamics of the tale is not adiabatic, that is, it is not completely controlled by the spectral energy cascade since the input of energy due to the nonlinear interaction competes with the energy input from wind (curve 2) and breaking dissipation (curves 3 and 4). The tail dissipation (curve 4) is small and concentrated in the vicinity of the cut wave number. The input and dissipation in the spectral tail are nearly in balance (curve 4) (Figure 8).

Figure 8.

The rates of transformation spectra multiplied by ${10}^9$ for the last period: (1) nonlinear interaction $N\left(r\right)$(1); (2) input energy; (3) breaking dissipation $D_b\left(r\right)$; (4) tail dissipations ${10}^{10}{D}_t\left(r\right)$; and (5) balance of all the terms.

6. Statistical properties of wave field

The phase-resolving modeling requires a higher computer capacity for calculations of any statistical characteristics of sea waves. In the course of simulations, 1.200 two-dimensional fields of the elevation and surface potential with the size 2048 × 1024 points were recorded. Following the solution of the 3-D equation for the velocity potential, these data allow us to reproduce any kinematic and dynamic characteristics of the three-dimensional structure of waves.

The most important statistical characteristics of wave field are mean $H_s$, variance V, skewness Sk, and kurtosis Ku calculated by the averaging over each of 1200 fields:

$$V=\overline{{\left(\eta -\overline{\eta }\right)}^2},\mathit{Sk}=\overline{{\left(\eta -\overline{\eta }\right)}^3}{V}^{-3/2},\mathit{Ku}==\overline{{\left(\eta -\overline{\eta }\right)}^4}{V}^{-2}-3.$$E38

The evolution of these characteristics in time is shown in Figure 9.

Figure 9.

Evolution of the statistical characteristics of elevation: (1) mean value $\overline{\eta }$, (2) variance $V$, (3) skewness Sk, and (4) kurtosis Ku.

The volume of the domain characterized by $\overline{\eta }$ is preserved with the accuracy of the order of ${10}^{-8}$. The variance $V$ is the potential energy that is growing up to the saturation. When the wave field is a superposition of a large number of linear waves, both the skewness and kurtosis are equal to zero. The skewness S characterizes asymmetry of the probability distribution indicating that the positive values of $\eta$ are larger than the negative ones, then $S>0$. The kurtosis $\mathit{Ku}$ is positive if the crests are sharper and the troughs are smoother than in the case of linear waves.

The probabilities of the geometrical characteristics (elevation, first and second derivatives over x) are shown in Figure 10. The elevation Z (normalized by the significant wave height) is characterized by asymmetry: the heights of waves are significantly larger than the depths of troughs, that is, the wave field is closer to the superposition of Stokes waves than to that of the harmonic modes. The distribution of slopes exhibits horizontal asymmetry: the negative slopes are larger than the positive ones, that is, the waves, on the average, are inclined in the direction of movements. The second derivative (curvilinearity) has the most striking tendency for asymmetry: the negative values corresponding to the sharpness of crests are much larger by absolute value than the positive values corresponding to the curvilinearity of troughs. This property of curvilinearity was used for the parameterizing of breaking. The limit value $Z_{\mathit{xx}}=-50$ was used as a criterion for the initiating of breaking (see Eq. (26)).

Figure 10.

Geometric characteristics of elevation: (1) probability of elevation $P\left(Z\right)$; (2) probability of slopes $P\left(Z_x\right)$; and (3) probability of curvilinearity $P\left(Z_{\mathit{xx}}\right)$.

The probability for three components of the surface velocity is given in Figure 11. The distributions of the vertical and transverse components of velocity are symmetrical. For a horizontal component, the values of positive fluctuations are considerably larger than the negative fluctuations. This effect cannot be explained by the influence of Stokes drift, which value for those specific conditions does not exceed ${10}^{-3}.$ The asymmetry of the probability distribution for the $u$-components is definitely connected with the asymmetry of the probability distribution for inclinations of surface (Figure 10, panel 2).

Figure 11.

Probability of the longitudinal $u$, transverse $v$, and vertical $w$ components of the surface velocity calculated for the last of nine periods corresponding approximately to the quasi-stationary regime.

The number of extreme waves with a high crest $Z_c/H_s>1.2$ is shown in Figure 12. Because such wave is not presented in each of the wave fields, the picture looks as discrete bars of different heights. The total number of values $Z/H_s>1.2$ is 17.214. The formally calculated probability of the values equals $0.67\cdot {10}^{-5}$. Note that the data on the probability of wave height contain uncertainty because it is not always clear which event should be considered as a single freak wave. The straightforward way suggests calculation of a portion of all the records including freak waves, out of the total volume of the data.

Figure 12.

The number of points where the nondimensional height $Z_c/H_s$ exceeds 1.2. The total number of values is $2.57\cdot {10}^9$.

However, some of the records can belong to the single moving freak waves. The cause of this uncertainty is the absence of a strict definition of freak wave being either a case or a process. The number of extreme waves grows with development of wave field.

The integral probability of the total wave height $Z_{\mathit{tc}}/H_s$, the wave height above mean level $Z_c/H_s$, and the depth of trough $Z_t/H_s$ are shown in Figure 13.

Figure 13.

The cumulative probability of crest-to-trough wave height $Z_{\mathit{tc}}/H_s$ (curve 3); crest height $Z_c/H_s$ (curve 2); and trough depth $-Z_t/H_s$ (curve 3). The number of points in each filed is equal to $2048\times 1024$. The number of fields is 1200.

Thin lines show that $Z_{\mathit{tc}}/H_s=2$ correspond approximately to $Z_c/H_s=1.2$ and $-Z_t/H_s=-.86$. It is worth to remind that here the nondimensional ‘extreme’ waves are considered. The true extreme waves are the product of the real wave field. The probability of real extreme waves can be estimated by multiplying the probability of the nondimensional wave by the probability of significant wave height.

The statistical connection between the total wave heights, crest heights, and trough depths is shown in Figure 14.

Figure 14.

Dependence of crest height $Z_c/H_s$ (top section) and depth trough $Z_t/H_s$(bottom section) on the total wave height $Z_{\mathit{tc}}/H_s$.

The dependences between these characteristics can be approximated by the formulas:

$$\begin{array}{l}{\tilde{Z}}_c=‐0.105+0.626{\tilde{\mathrm{Z}}}_{\mathit{tc}}+0.015{\tilde{\mathrm{Z}}}_{\mathit{tc}}^2\\ {\tilde{\mathrm{Z}}}_t=‐0.105‐0.374{\tilde{\mathrm{Z}}}_{\mathit{tc}}+0.015{\tilde{\mathrm{Z}}}_{\mathit{tc}}^2\end{array}$$E39

where the tilde denotes the normalizing by significant wave height $H_s$. Note that the first and third coefficients in (39) turned out to be a match. The correlation coefficient between ${\tilde{Z}}_t \text{and} {\tilde{Z}}_c$ is $-0.354$, while between ${\tilde{Z}}_t \text{and} {\tilde{\mathrm{Z}}}_{\mathit{tc}}$, it is −0.721, and between ${\tilde{Z}}_c \text{and} {\tilde{\mathrm{Z}}}_{\mathit{tc}}$, it is 0.903, that is, the correlation between the full wave height ${\tilde{Z}}_{\mathit{tc}}$ and the wave height above mean level ${\tilde{Z}}_c$ is so high that ${\tilde{Z}}_c$ can be used for identification of extreme waves.

The last characteristics that we consider here is the angle distribution of the spectral density. This characteristic can be described by the function $\mathrm{\Upsilon }\left(\omega /{\omega }_p\right)$ (see Ref. [60]).

$$\mathrm{\Upsilon }=\frac{\int S\left(\omega ,\psi \right)\left|\theta \right|\mathit{\text{dωdψ}}}{\int S\left(\omega ,\psi \right)\mathit{\text{dωdψ}}}$$E40

where the integrals are taken over the domain $\left\{\left(0<\omega <{\omega }_c\right),\left(-\pi /2<\psi <\pi /2\right)\right\}$. The value $\mathrm{\Upsilon }$ is weighted by the absolute spectrum value of wave direction. The wave spectra as the functions of frequency $\omega$ normalized by peak frequency ${\omega }_p$ for the first seven periods are shown in the upper panel of Figure 15.

Figure 15.

The shape of wave spectrum as a function of the nondimensional frequency $\omega /{\omega }_p$ (top panel) and a function $\mathrm{\Upsilon }$ (Eq. (40); ${\omega }_p$ is the frequency in the spectral peak.

The function $\mathrm{\Upsilon }\left(\omega /{\omega }_p\right)$ calculated for the same spectra is given in the bottom panel. As seen, the $\mathrm{\Upsilon }$ curves corresponding to different wave ages are close to each other. All of them have a sharp maximum at the frequencies below the spectral peak, a well-pronounced minimum in the spectral peak, and a relatively slow growth above the spectral peak. The decrease of $\mathrm{\Upsilon }$ at high frequencies is probably caused by the high-frequency dumping. The angle distribution was investigated in Refs. [56, 57, 58, 59, 60]. The approximations of $\mathrm{\Upsilon }\left(\omega /{\omega }_p\right)$ from the different sources collected in Ref. [61] show considerable scatter, but the general features are quite similar to those calculated in the current work. Note that the spectrum has undergone a long development; hence, the characteristics presented in Figure 14 were produced by the numerical model itself.

7. Conclusions

The paper is devoted to the wind wave simulations based on the initial equations of potential motion of fluid with a free surface. The system of equations includes the evolutionary kinematic and dynamic surface conditions and Laplace equation for the velocity potential. In this paper, a case of the double-periodic domain of infinite depth is considered. The construction of the exact numerical scheme for a long-term integration of these equations in the Cartesian coordinate system is impossible, since the surface moves between the grid knots. Instead, the system of the curvilinear coordinates (1) fitted with the surface is introduced. The main advantage of this coordinate system is that the surface coincides with a coordinate line $\zeta =0.$ The penalty follows immediately after turning the simple Cartesian coordinates into the curvilinear, nonstationary, and nonorthogonal coordinate system. Fortunately, the evolutionary Eqs. (4) and (5) become just slightly complicated, while Laplace equation transforms into the full elliptic equation. At each time step, these equations can be represented as Poisson equation with the right-hand side depending on the solution itself as well as on the metric coefficient. Since the norm of the right sight of the equation is usually small, the solution of Poisson equation can be found with the three-diagonal matrix algorithm and with iterations over the right-hand side. This procedure being formulated in the Fourier space is greatly simplified by the assumption of periodicity since in this case the derivatives over the horizontal coordinates are represented by the absolute value of wave number $\left|k\right|$ in the diagonal terms. When constructing a numerical scheme, we noticed that the significant simplification of the problem can be achieved by separation of the velocity potential into the linear and nonlinear components (see Ref. [7]). It is assumed that the linear component satisfies Laplace equations with the known solution. The equation for the nonlinear component can be obtained by extracting Laplace equation from the initial Poisson equation. Such procedure has a lot of advantages since the nonlinear component is on the average less by 1–2 decimal orders than the linear one. It means that for solution of the reduced Poisson equation the lesser number of levels in vertical, the lesser number of iterations and a smaller accuracy criterion can be used. The use of two components in the evolutionary equation does not seem to provide noticeable advantages; however, this way deserves further consideration.

The adiabatic version of the model was validated by simulation of a running Stokes wave with the steepness $\mathit{AK}=0.40$ in Ref. [7]. It was shown that the amplitudes of Stokes modes remain practically constant up to the accuracy of ${10}^{-7}$. The current version of the model after some technical improvements of the numerical scheme provides accuracy up to ${10}^{-12}$. Then, the adiabatic version of the model was used for reproduction of a quasi-stationary regime for investigation of the statistical properties of sea waves [1, 7, 8].

For calculations of development of wave field under the action of wind, it was necessary to include the algorithms for calculations of input and dissipation of energy. The scheme for calculation of the energy input was developed by Chalikov and Rainchik [3] on the basis of coupling the one-dimensional phase-resolving model and the two-dimensional boundary layer model with the second-order turbulence closure scheme. The parameterization suggested is still quasi-linear (similar to Miles’ scheme [30]), but in our opinion, it is the only scheme confirmed by the extended results of the numerical simulations. The theoretical and observational data on $\beta -$function are dramatically scattered (see Ref. [4], Figure 1).

For stabilization of the solution, the algorithm of high-frequency dumping in the Fourier space suggested by Chalikov and Sheinin [28] was used. The numerous attempts were made to improve that scheme (for example by reduction of the spectral interval of dumping) but without much success.

The most complicated problem is the parameterization of dissipation due to wave breaking. Such algorithm should not describe a process of breaking as it is, which within the frame of such model is impossible, but it should prevent the numerical instability that interrupts a run (see discussion in Ref. [26]). Currently, the algorithm used is very simple. It is based on the diffusion operator with a highly selective coefficient of ‘viscosity.’ It works satisfactorily, but we are far from thinking that it cannot be substantially improved or completely replaced by another one.

The results described in this paper show that the wave field development under the action of wind is reproduced quite realistically. The area of application of such models is very wide. Such modeling should be used for improvement of the algorithms of the energy input and dissipation. A model with the periodic boundary conditions can be used for the local interpretation of the spectral forecast in terms of real waves. The finite-difference version of the model can be used for simulation of wave regimes in the basins with real shapes and bathymetry (see, e.g., Ref. [5]).

Acknowledgments

Author would like to thank Mrs. O. Chalikova for her assistance in preparation of the manuscript. This research was performed in the framework of the state assignment of Russian Academy of Science (Theme No. 0149-2019-0015) supported in part 15 (Section 2) by RFBR (project No. 18-05-01122).