Numerical modeling of surface wave development under the action of wind

The numerical modeling of two-dimensional surface wave development under the action of wind is performed. The model is based on three-dimensional equations of potential motion with a free surface written in a surface-following nonorthogonal curvilinear coordinate system in which depth is counted from a moving surface. A three-dimensional Poisson equation for the velocity potential is solved iteratively. A Fourier transform method, a second-order accuracy approximation of vertical derivatives on a stretched vertical grid and fourth-order Runge–Kutta time stepping are used. Both the input energy to waves and dissipation of wave energy are calculated on the basis of earlier developed and validated algorithms. A one-processor version of the model for PC allows us to simulate an evolution of the wave field with thousands of degrees of freedom over thousands of wave periods. A long-time evolution of a two-dimensional wave structure is illustrated by the spectra of wave surface and the input and output of energy.

1 Introduction

The phase-resolving modeling of sea waves is the mathematical modeling of surface waves including explicit simulations of surface elevation and a velocity field evolution. As compared with spectral wave modeling, phase-resolving modeling is more general since it reproduces a real visible physical process and is based on well-formulated full equations. Phase-resolving models usually operate with a large number of degrees of freedom. In general, this method is more complicated and requires more computational resources. The simplest way to model like this is to calculate wave field evolution based on linear equations. Such an approach allows the reproduction of the main effects of the linear wave transformation due to the superposition of wave modes, reflections, refractions, etc. This approach is useful for many technical applications but it cannot reproduce a nonlinear nature of waves and the transformation of wave field due to the nonlinearity. Another example of a relatively simple object is a case of the shallow-water waves. The nonlinearity can be taken into account in the more sophisticated models derived from the fundamental fluid mechanics equations with some simplifications. The most popular approach is based on a nonlinear Schrödinger equation of different orders (see Dysthe, 1979) obtained by expansion of the surface wave displacement. This approach is also used for solving the problem of freak waves. The main advantage of a simplified approach is that it allows the reduction of a three-dimensional (3-D) problem to a two-dimensional one (or 2-D problem to 1-D problem). However, it is not always clear which of the nonrealistic effects is eliminated or included in the model after simplifications. This is why the most general approach being developed over the past years is based on the initial two-dimensional or three-dimensional equations (still potential). All the tasks based on these equations can be divided into two groups: the periodic and nonperiodic problems. An assumption of periodicity considerably simplifies construction of the numerical models, though such formulation can be applied to the cases when the condition of periodicity is acceptable, for example, when domain is considered to be a small part of a large uniform area. For the limited domains with no periodicity the problem becomes more complicated since the Fourier presentation cannot be used directly.

From the point of view of physics, the problem of phase-resolving modeling can be divided into two groups: the adiabatic and nonadiabatic modeling. A simple adiabatic model assumes that the process develops with no input or output of energy. Being not completely free of limitations, such a formulation allows the investigation of the wave motion on the basis of true initial equations. Including the effects of input energy and its dissipation is always connected with the assumptions that generally contradict the assumption of potentiality, i.e., the new terms added to the equations should be referred to as pure phenomenological. This is why the treatment of a nonadiabatic approach is often based on quite different constructions.

All of the phase-resolving models use the methods of computational mathematics and inherit all their advantages and disadvantages; i.e., on the one side, there is the possibility of a detailed description of the processes, and on the other side, there are a bunch of the specific problems connected with the computational stability, space and time resolution. The mathematical modeling produces tremendous volumes of information, the processing of which can be more complicated than the modeling itself.

The phase-resolving wave modeling takes a lot of computer time since it normally uses a surface-following coordinate system, which considerably complicates the equations. The most time-consuming part of the model is an elliptic equation for the velocity potential usually solved with iterations. Luckily, for a two-dimensional problem this trouble is completely eliminated by use of the conformal coordinates, reducing the problem to a one-dimensional system of equations which can be solved with high accuracy (Chalikov and Sheinin, 1998). For a three-dimensional problem, the reduction to a two-dimensional form is evidently impossible; hence, the solution of a 3-D elliptical equation for the velocity potential becomes an essential part of the entire problem. This equation is quite similar to the equation for pressure in a nonpotential problem. It follows that the 3-D Euler equations, being more complicated, can still be solved over the acceptable computer time.

There is a large volume of papers devoted to the numerical methods developed for the investigation of wave processes over the past decades. It includes a finite-difference method (Engsig-Karup et al., 2009, 2012), a finite-volume method (Causon et al., 2010), a finite-element method (Ma and Yan, 2010; Greaves, 2010), a boundary (integral) element method (Grue and Fructus, 2010), and spectral methods (Ducroset et al., 2007, 2012, 2016; Touboul and Kharif, 2010; Bonnefoy et al., 2010). These include a smoothed-particle hydrodynamics method (Dalrymple et al., 2010), a large-eddy simulation (LES) method (Issa et al., 2010; Lubin and Caltagirone, 2010), a moving particle semi-implicit method (Kim et al., 2014), a constrained interpolation profile method (Zhao, 2016), a method of fundamental solutions (Young et al., 2010) and a meshless local Petrov–Galerkin method (Ma and Yan, 2010). A fully nonlinear model should be applied to many problems. Most of the models were designed for engineering applications such as overturning waves, broken waves, waves generated by landslides, freak waves, solitary waves, tsunamis, violent sloshing waves, an interaction of extreme waves with beaches and an interaction of steep waves with the fixed structures or with different floating structures. The references given above make up less than 1 % of the publications on those topics.

A two-dimensional approach (like a conformal method) considers a strongly idealized wave field since even monochromatic waves in the presence of lateral disturbances quickly obtain a two-dimensional structure. The difficulty arising is not a direct result of the increase in the dimension. The fundamental complication is that the problem cannot be reduced to a two-dimensional problem, and even for the case of a double-periodic wave field, the problem of solution of a Laplace-like equation for the velocity potential arises. The majority of the models designed for investigation of the three-dimensional wave dynamics are based on simplified equations such as the second-order perturbation methods in which the higher-order terms are ignored. Overall, it is unclear which effects are missing in such simplified models.

The most sophisticated method is based on the full three-dimensional equations and surface integral formulations (Beale, 2001; Xue et al., 2001; Grilli et al., 2001; Clamond and Grue, 2001; Clamond et al., 2005, 2006; Fructus et al., 2005; Guyenne et al., 2006; Fochesato et al., 2006). A fully nonlinear model of three-dimensional water waves, which extends an approach suggested by Craig and Sulem (1993), was originally given in a two-dimensional setting. The model is based upon the Hamiltonian formulation (Zakharov, 1968), which allows the reduction of the problem of surface variable computation by introducing a Dirichlet–Neumann operator, which is expressed in terms of its Taylor series expansion in homogeneous powers of surface elevation. Each term in this Taylor series can be obtained from the recursion formula and efficiently computed using a fast Fourier transform.

The main advantage of the boundary integral equation methods (BIEMs) is that they are accurate and can describe highly nonlinear waves. A method of solution of the Laplace equation is based on the use of Green's function, which allows us to reduce a 3-D water wave problem to a 2-D boundary integral problem. The surface integral method is well suited for simulation of the wave effects connected with very large steepness, specifically, for investigation of the freak wave generation. These methods can be applied both to the periodic and nonperiodic flows. The methods do not impose any limitations on the wave steepness; thus they can be used for simulation of the waves that even approach breaking (Grilli et al., 2001) when the surface obtains a non-single value shape. The method allows us to take into account the bottom topography (Grue and Fructus, 2010) and investigate an interaction of waves with the fixed structures or with the freely responding floating structures (Liu et al., 2016; Gou et al., 2016).

However, the BIEM seems to be quite complicated and time consuming when applied to the long-term evolution of a multimode wave field in large domains. The simulation of the relatively simple wave fields illustrates an application of the method, and it is unlikely that the method can be applied to the simulation of the long-term evolution of a large-scale multimode wave field with a broad spectrum. An implementation of a multipole technique for a general problem of the sea wave simulation (Fochesato et al., 2006) can solve the problem but obviously leads to considerable algorithmic difficulties.

Currently, the most popular approach in the oceanography approach is a HOS (high-order scheme) model developed by Dommermuth and Yue (1987) and West et al. (1987). The HOS model is based on a paper by Zakharov (1968) in which a convenient form of the dynamic and kinematic surface conditions was suggested. The equations used by Zakharov were not intended for modeling, but rather for investigation of stability of the finite amplitude waves. In fact, a system of coordinates in which depth is counted from the surface was used, but the Laplace equation for the velocity potential was taken in its traditional form. However, the Zakharov's followers have accepted this idea literally. They used the two coordinate systems: a curvilinear surface-fitting system for the surface conditions and the Cartesian system for calculation of the surface vertical velocity. An analytical solution for the velocity potential in the Cartesian coordinate system is known. It is based on the Fourier coefficients on a fixed level, while the true variables are the Fourier coefficients for the potential on a free surface. Here a problem of transition from one coordinate system to another arises. This problem is solved by expansion of the surface potential into the Taylor series in the vicinity of the surface. The accuracy of this method depends on that of the representation of the exponential function with a finite number of the Taylor series. For the small-amplitude waves and for a narrow wave spectrum, such accuracy is evidently satisfactory. However, for the case of a broad wave spectrum that contains many wave modes, the order of the Taylor series should be high. The problem is now that the waves with high wave numbers are superposed over the surface of larger waves. Since the amplitudes of a surface potential attenuate exponentially, the amplitude of a small wave at a positive elevation increases, and conversely, it can approach zero at negative elevations. It is clear that such a setting of the HOS model cannot reproduce high-frequency waves, which actually reduces the nonlinearity of the model. This is why such a model can be integrated for long periods using no high-frequency smoothing. In addition, an accuracy of the calculation of a vertical velocity on the surface depends on full elevation at each point. Hence, the accuracy is not uniform along the wave profile. A substantial extension of the Taylor series can definitely result in numerical instability due to the occasional amplification of modes with high wave numbers. The authors of a surface integral method share a similar point of view (Clamond et al., 2005). We should note, however, that the comparison of the HOS method based on the West et al. (1987) approach using a method of the surface integral for an idealized wave field (Clamond et al., 2006) shows quite acceptable results. It was shown in the previous paper that a method suggested by Dommermuth et al. (1987) demonstrates poorer divergence of the expansion for the vertical velocity than the method by West et al. (1987). The HOS model has been widely used (for example, Tanaka, 2001; Toffoli et al., 2010; Touboul and Kharif, 2010) and it has shown its ability to efficiently simulate the wave evolution (propagation, nonlinear wave–wave interactions, etc.) in a large-scale domain (Ducrozet et al., 2007, 2012). It is obvious that the HOS model can be used for many practical purposes. Recently, Ecole Centrale Nantes, LHEEA Laboratory (CNRÑ) announced that the nonlinear wave models based on HOS are published as an open source (https://github.com/LHEEA/HOS-ocean/wiki, last access: 6 June 2018).

Opposite to the HOS method based on the analytical solution of the Laplace equation in Cartesian coordinates, a group of models is based on a direct solution of the equation for the velocity potential in the curvilinear coordinates (Engsig-Karup et al., 2009, 2012; Chalikov et al., 2014). The main advantage of a surface-following coordinate system is that a variable surface is mapped onto the fixed plane. Since the wave motion is very conservative, the highly accurate numerical schemes should be used for a good description of the nonlinearity and spectrum transformation. This most universal approach is being developed at the Technical University of Denmark (TUD) (see Engsig-Karup, 2009). Actually, the models ModelWave3D developed at TUD are targeted at the solution of a variety of problems, including such problems as the modeling of wave interaction with submerged objects as well as the simulation of wave regime in basins with a real shape and topography.

The model is based on the equations of a potential flow with a free surface. An effect of variable bathymetry is taken into account by using the so-called σ coordinate, (straightening out the bottom and surface). At vertical surfaces a normal derivative of the velocity potential is equal to zero. A flexible-order approximation for spatial derivatives is used. The most time-consuming part of this mode is a 3-D equation for the velocity potential. The strategy of the model development is directed at exploiting the architectural features of modern graphics processing units for the mixed precision computations. This approach is tested using a recently developed generic library for fast prototyping of PDE (partial differential equation) solvers. The new wave tool is applicable for solving and analyzing a variety of large-scale wave problems in coastal and offshore engineering. A description of the project and references can be found at the site https://www.researchgate.net/project/OceanWave3D-Open-Source-Marine-Hydrodynamics (last access: 6 June 2018).

A comparison of a ModelWave3D with a HOS model was presented by Ducrozet et al. (2012). It was shown that both models demonstrate high accuracy, while the HOS model shows a better performance. Note that the comparison of the speed of the models in this case is irrelevant since the ModelWave3D was designed for investigation of complicated processes, taking into account the real shape of a basin, variable depth and even the presence of engineering constructions. All these features are obviously not included in the HOS model.

The development of waves under the action of wind is a process that is difficult to simulate since surface waves are very conservative and change their energy for hundreds and thousands of periods. This is why the most popular method is spectral modeling. Waves as physical objects in this approach are actually absent since an evolution of the spectral distribution of wave energy is simulated. The description of input and dissipation in this approach is not directly connected with the formulation of the problem, but rather it is adopted from other branches of wave theory in which waves are the objects of investigation. However, the spectral approach was found to be the only method capable of describing the space and time evolution of wave field in the ocean. The phase-resolving models (or “direct” models) designed for reproducing the waves themselves cannot compete with the spectral models since the typical size of the domain in such models does not exceed several kilometers. Such a domain includes just several thousands of large waves. Nevertheless, the direct wave modeling plays an ever-increasing role in geophysical fluid dynamics because it gives the possibility of investigating the processes which cannot be reproduced with spectral models. One such problem is that of an extreme wave generation (Chalikov, 2009; Chalikov and Babanin, 2016a). Direct modeling is also a perfect instrument for the development of parameterization of physical processes for spectral wave models. In addition, such models can be used for direct simulation of wave regimes of small water basins, for example, port harbors. Other approaches of direct modeling are discussed in Chalikov et al. (2014) and Chalikov (2016).

Until recently the direct modeling was used for reproduction of a quasi-stationary wave regime when the wave spectrum did not change significantly. A unique example of the direct numerical modeling of a surface wave evolution is given in Chalikov and Babanin (2014), in which the development of a wave field was calculated with the use of a two-dimensional model based on the full potential equations written in the conformal coordinates. The 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 an evolution of wave spectrum under the action of wind. However, the strictly one-dimensional (unidirected) waves are not realistic; hence, a full problem of wave evolution should be formulated on the basis of the three-dimensional equations. An example of such modeling is given in the current paper.

2 Equations

Let us introduce a nonstationary surface-following nonorthogonal coordinate system:

$$\begin{array}{}\text{(1)}& \begin{array}{cccc}\mathit{\xi }=x,& \mathit{\vartheta }=y,& \mathit{\zeta }=z-\mathit{\eta }(\mathit{\xi },\mathit{\vartheta },\mathit{\tau }),& \mathit{\tau }=t\end{array},\end{array}$$

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

$$\begin{array}{}\text{(2)}& \mathit{\eta }\left(\mathit{\xi },\mathit{\vartheta },\mathit{\tau }\right)=\sum _{-M_x<k<M_x}\sum _{-M_y<l<M_y}{h}_{k,l}\left(\mathit{\tau }\right){\mathrm{\Theta }}_{k,l},\end{array}$$

where k and l are the components of a wave number vector k, h_k,l(τ) are Fourier amplitudes for elevations $\mathit{\eta }\left(\mathit{\xi },\mathit{\vartheta },\mathit{\tau }\right)$, M_x and M_y are the numbers of modes in the directions ξ and ϑ, respectively, and Θ_k,l are the Fourier expansion basis functions, represented as the matrix

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

The 3-D equations of potential waves in the system of coordinates (1) at ζ≤0 take the following form:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathit{\tau }}=-{\mathit{\eta }}_{\mathit{\xi }}{\mathit{\phi }}_{\mathit{\xi }}-{\mathit{\eta }}_{\mathit{\vartheta }}{\mathit{\phi }}_{\mathit{\vartheta }}+\left(\mathrm{1}+{\mathit{\eta }}_{\mathit{\xi }}^{\mathrm{2}}+{\mathit{\eta }}_{\mathit{\vartheta }}^{\mathrm{2}}\right){\mathrm{\Phi }}_{\mathit{\varsigma }},\text{(5)}& {\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathit{\tau }}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\phi }}_{\mathit{\xi }}^{\mathrm{2}}+{\mathit{\phi }}_{\mathit{\vartheta }}^{\mathrm{2}}-\left(\mathrm{1}+{\mathit{\eta }}_{\mathit{\xi }}^{\mathrm{2}}+{\mathit{\eta }}_{\mathit{\vartheta }}^{\mathrm{2}}\right){\mathrm{\Phi }}_{\mathit{\zeta }}^{\mathrm{2}}\right)-\mathit{\eta }-p,\text{(6)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathit{\xi }\mathit{\xi }}+{\mathrm{\Phi }}_{\mathit{\vartheta }\mathit{\vartheta }}+{\mathrm{\Phi }}_{\mathit{\zeta }\mathit{\zeta }}=\mathrm{{\rm Y}}\left(\mathrm{\Phi }\right),\end{array}$$

where Υ is the operator:

$$\begin{array}{}\text{(7)}& \mathrm{{\rm Y}}\left(\right)=\mathrm{2}{\mathit{\eta }}_{\mathit{\xi }}({)}_{\mathit{\xi }\mathit{\zeta }}+\mathrm{2}{\mathit{\eta }}_{\mathit{\vartheta }}({)}_{\mathit{\vartheta }\mathit{\zeta }}+\left({\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}+{\mathit{\eta }}_{\mathit{\vartheta }\mathit{\vartheta }}\right)({)}_{\mathit{\zeta }}-\left({\mathit{\eta }}_{\mathit{\xi }}^{\mathrm{2}}+{\mathit{\eta }}_{\mathit{\vartheta }}^{\mathrm{2}}\right)({)}_{\mathit{\zeta }\mathit{\zeta }}.\end{array}$$

Capital fonts Φ are used for domain ζ<0 while the lower case φ refers to ζ=0. A term p in Eq. (5) described the pressure on the surface ζ=0.

It is suggested in Chalikov et al. (2014) that it is convenient to represent the velocity potential φ as a sum of two components, i.e., an analytical (linear) component $\overline{\mathrm{\Phi }}\left(\overline{\mathit{\phi }}=\overline{\mathrm{\Phi }}\left(\mathit{\xi },\mathit{\vartheta },\mathrm{0}\right)\right)$ and an arbitrary (nonlinear) component $\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\stackrel{\mathrm{̃}}{\mathit{\phi }}=\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\left(\mathit{\xi },\mathit{\vartheta },\mathrm{0}\right)\right)$:

$$\begin{array}{}\text{(8)}& \mathit{\phi }=\stackrel{\mathrm{‾}}{\mathit{\phi }}+\stackrel{\mathrm{̃}}{\mathit{\phi }},\mathrm{\Phi }=\stackrel{\mathrm{‾}}{\mathrm{\Phi }}+\stackrel{\mathrm{̃}}{\mathrm{\Phi }}.\end{array}$$

The analytical component $\overline{\mathrm{\Phi }}$ satisfies the Laplace equation

$$\begin{array}{}\text{(9)}& {\overline{\mathrm{\Phi }}}_{\mathit{\xi }\mathit{\xi }}+{\overline{\mathrm{\Phi }}}_{\mathit{\vartheta }\mathit{\vartheta }}+{\overline{\mathrm{\Phi }}}_{\mathit{\zeta }\mathit{\zeta }}=\mathrm{0},\end{array}$$

with a known solution:

$$\begin{array}{}\text{(10)}& \overline{\mathrm{\Phi }}(\mathit{\xi },\mathit{\vartheta },\mathit{\zeta },\mathit{\tau })=\sum _{k,l}{\overline{\mathit{\phi }}}_{k,l}\left(\mathit{\tau }\right)\exp \left(\left|k\right|\mathit{\zeta }\right){\mathrm{\Theta }}_{k,l}\end{array}$$

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

$$\begin{array}{}\text{(11)}& \begin{array}{l}\mathit{\varsigma }=\mathrm{0}:\overline{\mathrm{\Phi }}=\overline{\mathit{\phi }}\\ \mathit{\varsigma }\to -\mathrm{\infty }:{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_{\mathit{\zeta }}\to \mathrm{0}\end{array}.\end{array}$$

The nonlinear component satisfies an equation:

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_{\mathit{\xi }\mathit{\xi }}+{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_{\mathit{\vartheta }\mathit{\vartheta }}+{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_{\mathit{\zeta }\mathit{\zeta }}=\mathrm{{\rm Y}}\left(\stackrel{\mathrm{̃}}{\mathrm{\Phi }}\right)+\mathrm{{\rm Y}}\left(\overline{\mathrm{\Phi }}\right).\end{array}$$

Equation (12) is solved with the boundary conditions

$$\begin{array}{}\text{(13)}& \begin{array}{l}\mathit{\varsigma }=\mathrm{0}:\stackrel{\mathrm{̃}}{\mathrm{\Phi }}=\mathrm{0}\\ \mathit{\varsigma }\to -\mathrm{\infty }:{\stackrel{\mathrm{̃}}{\mathrm{\Phi }}}_{\mathit{\zeta }}\to \mathrm{0}\end{array}.\end{array}$$

The derivatives of a linear component $\overline{\mathrm{\Phi }}$ in Eq. (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 a stretched staggered grid defined by the relation $\mathrm{\Delta }{\mathit{\zeta }}_{\mathrm{j}+\mathrm{1}}=\mathit{\chi }\mathrm{\Delta }{\mathit{\zeta }}_{\mathrm{j}}$ (Δζ is a vertical step, while j = 1 at the surface). A stretched grid provides an increase in accuracy of approximation for the exponentially decaying modes. The values of a stretching coefficient χ lie within the interval 1.01–1.20. A finite-difference second-order approximation of the vertical operators in Eq. (12) on a nonuniform vertical grid is quite straightforward. Equation (12) is solved as the Poisson equations with the iterations over the right-hand side. At each time step, the iterations start with the right-hand side calculated at the previous time step. The initial elevation was generated as a superposition of linear waves corresponding to a JONSWAP spectrum (Hasselmann et al., 1973) with random phases. The initial Fourier amplitudes for the surface potential were calculated using the formulas of the linear wave theory. A detailed description of the scheme and its validation is given in Chalikov et al. (2014) and Chalikov (2016).

Equations (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^{\mathrm{1}/\mathrm{2}}{g}^{-\mathrm{1}/\mathrm{2}}$ and the velocity potential $L^{\mathrm{3}/\mathrm{2}}{g}^{\mathrm{1}/\mathrm{2}}$ (g is an acceleration of gravity). The pressure is normalized by water density, so that the pressure scale is Lg. Equations (4)–(6) are self-similar to the transformation with respect to L. The dimensional size of the domain is 2πL, so the scaled size is 2π. 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π and a factor $r=L_x/L_y$ is introduced into the definition of a differential operator in the Fourier space.

3 Energy input and dissipation

The energy input to waves is described by a pressure term p in a dynamic boundary condition (Eq. 5). The tangent stress on the surface cannot be taken into account in the potential formulation. The dissipation cannot be described either with use of the potential equations, but for a realistic description of wave dynamics, the dissipation of wave energy should be taken into account, i.e., we should include additional terms in Eqs. (4) and (5), which, strictly speaking, contradict an assumption of potentiality.

3.1 Energy input from wind

According to the linear theory (Miles, 1957), the Fourier components of surface pressure p are connected with those of the surface elevation through the following expression:

$$\begin{array}{}\text{(14)}& p_{k,l}+i{p}_{-k,-l}={\displaystyle \frac{{\mathit{\rho }}_a}{{\mathit{\rho }}_w}}\left({\mathit{\beta }}_{k,l}+i{\mathit{\beta }}_{-k,-l}\right)\left(h_{k,l}+i{h}_{-k,-l}\right),\end{array}$$

where $h_{k,l},h_{-k,-l},{\mathit{\beta }}_{k,l},{\mathit{\beta }}_{-k,-l}$ are real and imaginary parts of elevation η and the so-called β function (i.e., Fourier coefficients at COS and SIN, respectively); ρ_a∕ρ_w is a ratio of the air and water densities. Equation (14) is a standard presentation of pressure above a multimode surface. It means that every wave mode with amplitude ${\left(h_{k,l}^{\mathrm{2}}+h_{-k,-l}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}$ initiates the pressure mode with amplitude ${\left(p_{k,l}^{\mathrm{2}}+p_{-k,-l}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}}$ shifted off the phase of a wave mode by angle $\mathit{\alpha }=\mathrm{atan}\frac{{\mathit{\beta }}_{-k,-l}}{{\mathit{\beta }}_{k,l}}$. Both coefficients in Eq. (14) are a function of the ratio of wind velocity at half the mode of the length height ${\mathit{\lambda }}_{k,l}/\mathrm{2}$ to the virtual phase velocity. Hence, for derivation of the shape of the β function it is necessary to simultaneously measure the wave surface elevation and the nonstatic pressure on the surface. Measurement of surface pressure is a very difficult problem since measurements should be carried out very close to a moving surface, preferably, with a surface-following sensor. Such measurements are performed quite seldom, especially, in the field. The measurements were carried out for the first time by a team of authors in both laboratory and field (Snyder et al., 1981; Hsiao and Shemdin, 1983; Hasselmann and Bösenberg, 1991; Donelan et al., 2005, 2006). The data obtained in this way allowed the construction of an imaginary part of the β function used in some versions of wave forecasting models (Rogers et al., 2012). Such measurements and their processing are quite complicated since the wave-produced pressure fluctuations are masked by the turbulent pressure fluctuations. The second method of the β function evaluation is based on the results of numerical investigations of the statistical structure of a 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 using numerical models, not experimentally (Gent and Taylor, 1976; Riley et al., 1982; Al-Zanaidi and Hui, 1984). This method was being developed beginning from Chalikov (1978, 1986), followed by Chalikov and Makin (1991), Chalikov and Belevich (1992) and Chalikov (1995). The results were implemented in a WAVEWATCH model, i.e., a third-generation wave forecast model (Tolman and Chalikov, 1996), and thoroughly validated against the experimental data in the course of developing WAVEWATCH III (Tolman et al., 2014). This method was later improved on the basis of a more advanced coupled modeling of waves and boundary layer (Chalikov and Rainchik, 2010), while the β function used in WAVEWATCH III was corrected and extended up to the high frequencies. A direct calculation of the energy input to waves requires both the real and imaginary parts of the β function. The total energy input to waves depends on an imaginary part of the β function, while the moments of a higher order depend on both the imaginary and real parts of β. This is why a full approximation constructed in Chalikov and Rainchik (2010) was used in the current work. Note that in a range of relatively low frequencies the new method is very close to the scheme implemented in WAVEWATCH III.

It is a traditional suggestion that both coefficients are the functions of the virtual nondimensional frequency $\mathrm{\Omega }={\mathit{\omega }}_{k,l}U\cos \mathit{\psi }=U/c_{k,l}\cos \mathit{\psi }$ (where ω_k,l and U are the nondimensional radian frequency and wind speed, respectively; c_k,l is a phase speed of the kth mode; ψ is an angle between the wind and wave mode directions). Most of the schemes for calculations of the β function consider a relatively narrow interval of the nondimensional frequencies Ω. In the current work, the range of frequencies cover an interval $\left(\mathrm{0}<\mathrm{\Omega }<\mathrm{10}\right)$, and occasionally the values of Ω>10 can appear. This is another reason why the function derived in Chalikov and Rainchik (2010) through the coupled simulations of waves and the boundary layer is used here. The wave model is based on the potential equations for a flow with a free surface, extended with an algorithm for the breaking dissipation (see below a description of the breaking dissipation parameterization). The wave boundary layer (WBL) model is based on Reynolds equations closed with a K−ε scheme; the solutions for air and water are matched through the interface. The β function was used for the evaluation of accuracy of the surface pressure p calculations. The shape of the β function connecting surface elevations and surface pressure is studied up to the high nondimensional wave frequencies in both positive and negative (i.e., for wind opposite to waves) domains. The data on the β function exhibit wide scatter, but since the volume of the data was quite large (47 long-term numerical runs allowed us to generate about 1 400 000 values of β), the shape of the β function was defined with satisfactory accuracy up to the very high nondimensional frequencies $\left(-\mathrm{50}<\mathrm{\Omega }<\mathrm{50}\right)$. As a result, the data on the β-function in such a broad range allow us to calculate the wave drag up to very high frequencies and to explicitly divide the fluxes of energy and momentum transferred by the pressure and molecular viscosity. This method is free of arbitrary assumptions on a drag coefficient C_d, and, conversely, such calculations allow the investigation of the nature of the wave drag (see Ting et al., 2012).

The most reliable data on β function are concentrated in the interval $-\mathrm{10}<\mathrm{\Omega }<\mathrm{10}$ (negative values of Ω correspond to the wave modes running against wind). The real and imaginary parts of the β function are shown in Fig. 1. It is a corrected version of an approximation given in Chalikov and Rainchik (2010) in which the data at negative Ω were interpreted erroneously. In the current calculations the modes running against wind are absent. The function β can be approximated by the formulas

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{k,l}=\left\{\begin{array}{cc}{\mathit{\beta }}_{\mathrm{0}}+a_{\mathrm{0}}\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{0}}\right)+a_{\mathrm{1}}{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{0}}\right)}^{\mathrm{2}}& {\mathrm{\Omega }}_{\mathrm{0}}<\mathrm{\Omega }\\ {\mathit{\beta }}_{\mathrm{0}}+a_{\mathrm{0}}\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{0}}\right)-a_{\mathrm{1}}{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{0}}\right)}^{\mathrm{2}}& \mathrm{\Omega }<{\mathrm{\Omega }}_{\mathrm{0}}\end{array}\right.,\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{-k,-l}=\left\{\begin{array}{cc}{\mathit{\beta }}_{\mathrm{1}}+a_{\mathrm{3}}\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{2}}\right)& \mathrm{\Omega }<{\mathrm{\Omega }}_{\mathrm{2}}\\ a_{\mathrm{2}}{\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{1}}\right)}^{\mathrm{2}}& {\mathrm{\Omega }}_{\mathrm{2}}<\mathrm{\Omega }<{\mathrm{\Omega }}_{\mathrm{3}}\\ {\mathit{\beta }}_{\mathrm{1}}-a_{\mathrm{3}}\left(\mathrm{\Omega }-{\mathrm{\Omega }}_{\mathrm{3}}\right)& {\mathrm{\Omega }}_{\mathrm{3}}<\mathrm{\Omega }\end{array}\right.,\end{array}$$

where the coefficients are $a_{\mathrm{1}}=\mathrm{0.09476},a_{\mathrm{2}}=-\mathrm{0.3718},a_{\mathrm{3}}=\mathrm{14.80},{\mathit{\beta }}_{\mathrm{0}}=-\mathrm{0.02},\mathit{\beta }=-\mathrm{148}.a_{\mathrm{0}}=\mathrm{0.02277},{\mathrm{\Omega }}_{\mathrm{0}}=\mathrm{0.02277},{\mathrm{\Omega }}_{\mathrm{1}}=\mathrm{1.20},{\mathrm{\Omega }}_{\mathrm{2}}=-\mathrm{18.8},{\mathrm{\Omega }}_{\mathrm{3}}=\mathrm{21.2}$.

It was indicated above that an initial wave field is assigned as a superposition of linear modes whose amplitudes are calculated with a JONSWAP spectrum with an initial peak wave number $k_{\mathrm{p}}^{\mathrm{0}}=\mathrm{100}$. An initial value $U/c_{\mathrm{p}}^{\mathrm{0}}=\mathrm{6}$ was chosen, i.e., a ratio of the nondimensional wind speed at the height of one-half the initial peak wave length ${\mathit{\lambda }}_{\mathrm{0}}/\mathrm{2}=\mathrm{2}\mathit{\pi }/\mathrm{100}$, and the phase speed $c_{\mathrm{p}}^{\mathrm{0}}={\left(k_{\mathrm{p}}^{\mathrm{0}}\right)}^{-\mathrm{1}/\mathrm{2}}$ is equal to 6. Such a high ratio corresponds to the initial stages of wave development. The wind velocity $\mathrm{6}{c}_{\mathrm{p}}^{\mathrm{0}}$ remains constant throughout the integration. The values of Ω for other wave numbers are calculated by assuming that the wind profile is logarithmic:

$$\begin{array}{}\text{(17)}& {\mathrm{\Omega }}_{k,l}={\displaystyle \frac{U}{c_{k,l}}}\ln {\displaystyle \frac{{\mathit{\lambda }}_{k,l}}{\mathrm{2}{z}_{\mathrm{0}}}}{\left(\ln {\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{0}}}{\mathrm{2}{z}_{\mathrm{00}}}}\right)}^{-\mathrm{1}}\cos {\mathit{\psi }}_{k,l},\end{array}$$

where z₀₀ is an effective nondimensional roughness for the initial wind profile, while z₀ is an actual roughness parameter that depends on the energy in a high-frequency part of the spectrum and on the wind profile. We call it “effective” since very close to the surface the wind profile is not logarithmic (Chalikov, 1995; Tolman and Chalikov, 1996; Chalikov and Rainchik, 2010). The value of this parameter depends on the wind velocity and energy in a high-wave-number interval of the wave spectrum, as well as on the length scale of the problem. All these effects are possible to include by matching the wave model with a one-dimensional WBL model (Ting et al., 2012). Here, a simplified scheme for the roughness parameter is chosen. It is well known that the roughness parameter (as well as a drag coefficient) increases with a decrease in the inverse wave age. In our case the wind speed is fixed, and dependence for a nondimensional roughness parameter is constructed on the basis of the results obtained in Chalikov and Rainchik (2010):

$$\begin{array}{}\text{(18)}& z_{\mathrm{0}}=\mathrm{15}{z}_{\mathrm{00}}\mathrm{\Omega },\end{array}$$

where $z_{\mathrm{00}}={\mathrm{10}}^{-\mathrm{3}}$ is an initial value of the roughness parameter. Equation (18) approximates the dependence of the effective roughness at the stage of wave development. Note that the results are not sensitive to the variation of the roughness parameter within reasonable limits.

Figure 1 Real (dashed curve) and imaginary (solid curve) parts of the β function.

3.2 High-wave-number energy dissipation

A nonlinear flux of energy directed to the small wave numbers produces the downshifting of the spectrum, while an opposite flux forms the shape of the spectral tail. The second process can produce the accumulation of energy near a “cut” wave number. Both processes become more intensive with an increase in the energy input. The growth of amplitudes at high wave numbers is followed by that of the local steepness and numerical instability. This well-known phenomenon in the numerical fluid mechanics is eliminated by use of a highly selective filter simulating the nonlinear viscosity. To support stability, additional terms are included in the right-hand sides of Eqs. (4) and (5):

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\eta }}_{k,l}}{\partial \mathit{\tau }}}=E_{k,l}-{\mathit{\mu }}_{k,l}{\mathit{\eta }}_{k,l},\text{(20)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\phi }}_{k,l}}{\partial \mathit{\tau }}}=F_{k,l}-{\mathit{\mu }}_{k,l}{\mathit{\phi }}_{k,l}.\end{array}$$

E_k,l and F_k,l are the Fourier amplitudes of the right-hand sides of Eqs. (4) and (5) while a factor μ_k,l is calculated using the formula

$$\begin{array}{}\text{(21)}& {\mathrm{\mu }}_{k,l}=\left\{\begin{array}{cc}\mathrm{0}& \left|k\right|<k_d\\ c_m{k}_{\mathrm{0}}{\left({\displaystyle \frac{\left|k\right|-k_d}{\left(k_{\mathrm{0}}-k_d\right)}}\right)}^{\mathrm{2}}& k_d\le \left|k\right|\le k_{\mathrm{0}}\\ c_m{k}_{\mathrm{0}}& \left|k\right|>k_{\mathrm{0}}\end{array}\right.,\end{array}$$

where k and l are components of wave number |k|, while the coefficients k_d and k₀ are defined by the expression:

$$\begin{array}{}\text{(22)}& k_d=d_m^{\mathrm{2}}{M}_x{M}_y{\left({\left(l{\left|k\right|}^{-\mathrm{1}}{d}_m{M}_x\right)}^{\mathrm{2}}+{\left(k{\left|k\right|}^{-\mathrm{1}}{d}_m{M}_y\right)}^{\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(23)}& k_{\mathrm{0}}=M_x{M}_y{\left({\left(l{\left|k\right|}^{-\mathrm{1}}{M}_x\right)}^{\mathrm{2}}+{\left(k{\left|k\right|}^{-\mathrm{1}}{M}_y\right)}^{\mathrm{2}}\right)}^{-\mathrm{1}/\mathrm{2}},\end{array}$$

where c_m=0.1 and d_m=0.75. The expressions (19)–(21) can be interpreted in a straightforward way: the value of µ_k,l is equal to zero inside the ellipse with semiaxes d_mM_x and d_mM_y; then it grows linearly with |k| up to the value c_m and is equal to c_m outside the outer ellipse. This method of filtration that we call “tail dissipation” was developed and validated with a conformal model by Chalikov and Sheinin (1998). The sensitivity of the results to the parameters in Eqs. (21)–(23) is not high. The aim of the algorithm is to support smoothness and monotonicity of the wave spectrum within a high wave number range. Since the algorithm affects the amplitudes of small modes, it actually does not reduce the total energy, though it efficiently prevents development of the numerical instability. Note that no long-term calculations can be performed without tail dissipation eliminating the development of numerical instability at high wave numbers.

3.3 Dissipation due to wave breaking

The main process of wave dissipation is wave breaking. This process is taken into account in all the spectral wave forecasting models similar to WAVEWATCH (see Tolman and Chalikov, 1996). 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 spectral models in a distorted form. A real breaking occurs in relatively narrow areas of the physical space; however, a spectral image of such breaking is stretched over the entire wave spectrum, while in reality the breaking decreases height and energy of dominant waves. This contradiction occurs because the waves in spectral models are assumed to be the linear ones, while in fact the breaking occurs in the physical space with a nonlinear sharp wave, usually composed of several modes. However, progress has been gradually made in spectral wave modeling over the past decade. It became clear that state-of-the-art wave models should account for the threshold behavior of the dominant wave breaking, i.e., waves will not break unless their steepness exceeds the threshold (Alves and Banner, 2003; Babanin et al., 2010).

The mechanics of wave breaking at a developed wave spectrum differs from that in a wave field represented by the few modes normally considered in many theoretical and laboratory investigations (e.g., Alberello et al., 2018). Since the breaking in laboratory conditions is initiated by special assignment of amplitudes and phases, it cannot be similar to the breaking in natural conditions. To some degree the wave breaking is similar to the development of an extreme wave that appears suddenly with no pronounced prehistory (Chalikov and Babanin, 2016a, b). There are no signs of modulational instability in both phenomena, which suggests a process of energy consumption from other modes. The evolution leading to the breaking or “freaking” seems just opposite: the full energy of a main wave remains nearly constant while the columnar energy is focused around the crest of this wave, which becomes sharper and unstable. Probably even more frequent cases of wave breaking and extreme wave appearance can be explained by a local superposition of several modes.

The instability of interface leading to the breaking is an important and poorly developed problem of fluid mechanics. In general, this essentially nonlinear process should be investigated for a two-phase flow. Such an approach was demonstrated, for example, by Iafrati (2009). However, progress in solving this highly complicated problem is slow.

A problem of the breaking parameterization includes two points: (1) establishing of a criterion of the breaking onset and (2) development of an algorithm of the breaking parameterization. The problem of breaking is discussed in detail in Babanin (2011). Chalikov and Babanin (2012) performed a numerical investigation of the processes leading to the breaking. It was found that a clear predictor of the breaking formulated in dynamical and geometrical terms probably does not exist. The most evident criterion of the breaking is the breaking itself, i.e., the process when some part of the upper portion of a sharp wave crest is falling down. This process is usually followed by separation of the detached volume of liquid into the water and air phases. Unfortunately, there is no possibility of describing this process within the scope of the potential theory.

Some investigators suggest using a physical velocity approaching the rate of surface movement in the same direction as a criterion of the breaking onset. This is incorrect since a kinematic boundary condition suggests that these quantities are exactly equal to each other. It is quite clear that the onset of breaking can be characterized by the appearance of a non-single-value piece of surface. This stage can be investigated with a two-dimensional model, which due to a high flexibility of the conformal coordinates allows us to reproduce a surface with an inclination in the Cartesian coordinates exceeding 90 degrees. (In the conformal coordinates the dependence of elevation on a curvilinear coordinate is always a single value). The duration of this stage is extremely short, the calculations being always interrupted by the numerical instability with sharp violation of the conservation laws (constant integral invariants, i.e., full energy and volume) and strong distortion of the local structure of flow. The numerous numerical experiments with a conformal model showed that after the appearance of a non-single value the model never returns to stability. However, the introduction of a non-single surface as a criterion of the breaking instability even in a conformal model is impossible since a behavior of the model at a critical point is unpredictable, and the run is most likely to be terminated, no matter what kind of parameterization of breaking is introduced. It means that even in a precise conformal model the stabilization of the solution should be initiated prior to the breaking.

A consideration of an exact criterion for the breaking onset for the models using transformation of the coordinate type of point (1) is useless since the numerical instability in such models arises not because of the approach of breaking but because of the appearance of large local steepness. The multiple experiments with a direct 3-D wave model show that the appearance of the local steepness $max(\partial \mathit{\eta }/\partial x,\partial \mathit{\eta }/\partial y)$ exceeding ≈2 (that corresponds to a slope of about 60 degrees) is always followed by the numerical instability but the instability can happen well before reaching this value. The decrease in a time step does not make any effect. As seen, a surface with such a slope is very far from being a vertical wall, when the real breaking starts. However, an algorithm for the breaking parameterization must prevent numerical instability. The situation is similar to the numerical modeling of turbulence (LES technique) in which a local highly selective viscosity is used to prevent the appearance of too large local gradients of the velocity. A description of the breaking in the direct wave modeling should satisfy the following conditions. (1) It should prevent the onset of instability at each point of half a million grid points over more than 100 thousand time steps. (2) It should describe in a more or less realistic way the loss of kinetic and potential energies with preservation of balance between them. (3) It should preserve the volume. It was suggested in Chalikov (2005) 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:

$$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle }{\mathit{\eta }}_{\mathit{\tau }}=E_{\mathit{\eta }}+J^{-\mathrm{1}}\left({\displaystyle \frac{\partial }{\partial \mathit{\xi }}}{B}_{\mathit{\xi }}{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\xi }}}+{\displaystyle \frac{\partial }{\partial \mathit{\vartheta }}}{B}_{\mathit{\vartheta }}{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\vartheta }}}\right),\text{(25)}& {\displaystyle }& {\displaystyle }{\mathit{\phi }}_{\mathit{\tau }}=F_{\mathit{\phi }}+J^{-\mathrm{1}}\left({\displaystyle \frac{\partial }{\partial \mathit{\xi }}}{B}_{\mathit{\xi }}{\displaystyle \frac{\partial \mathit{\phi }}{\partial \mathit{\xi }}}+{\displaystyle \frac{\partial }{\partial \mathit{\vartheta }}}{B}_{\mathit{\vartheta }}{\displaystyle \frac{\partial \mathit{\phi }}{\partial \mathit{\vartheta }}}\right),\end{array}$$

where F_η and F_φ are the right-hand sides of Eqs. (4) and (5) including the terms introduced in terms of Fourier coefficients by Eqs. (19)–(23); B_ξ and B_ϑ are diffusion coefficients. It was suggested in the first versions of the scheme that a diffusion coefficient depends on a local slope; however, such a scheme did not prove to be very reliable since it did not prevent all of the events of the numerical instability. A scheme based on the calculation of the local curvilinearity η_ξξ and η_ϑϑ turned out to be a lot more robust. The calculations of 75 different runs were performed with a full 3-D model in Chalikov et al. (2014) over the period of t=350 (70 000 time steps). The total number of values used for the calculations of dependence in Fig. 2 (thick curve) is about 6 billion. The normal probability calculated with the same dispersion is shown by a thin curve.

It is seen that the probability of large negative values of the curvilinearity is calculated over an ensemble of linear modes by orders larger than the probability with the spectra generated by a 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 breaking approach. Coefficients B_ξ and B_ϑ depend nonlinearly on the curvilinearity:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }{B}_{\mathit{\xi }}=\left\{\begin{array}{cc}\mathrm{\Delta }\mathit{\xi }{C}_B{\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{\mathrm{2}}& {\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}<{\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{cr}\\ \mathrm{0}& {\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}\ge {\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{cr}\end{array}\right.,\text{(27)}& {\displaystyle }& {\displaystyle }{B}_{\mathit{\vartheta }}=\left\{\begin{array}{cc}\mathrm{\Delta }\mathit{\vartheta }{C}_B{\mathit{\eta }}_{\mathit{\vartheta }\mathit{\vartheta }}^{\mathrm{2}}& {\mathit{\eta }}_{\mathit{\vartheta }\mathit{\vartheta }}<{\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{cr}\\ \mathrm{0}& {\mathit{\eta }}_{\mathit{\vartheta }\mathit{\vartheta }}\ge {\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{cr}\end{array}\right.,\end{array}$$

where Δξ and Δζ are the horizontal steps in the x and y directions in a grid space, and the coefficients are C_B=2.0 and ${\mathit{\eta }}_{\mathit{\xi }\mathit{\xi }}^{cr}={\mathit{\eta }}_{\mathit{\vartheta }\mathit{\vartheta }}^{cr}=-\mathrm{50}$. Equations (24)–(27) do not change the volume and decrease the local potential and kinetic energy. It is assumed that the lost momentum and energy are transferred to the current and turbulence (see Chalikov and Belevich, 1992). In addition, the energy also goes to other wave modes. The choice of parameters in Eqs. (24)–(27) is based on simple considerations: a local piece of surface can closely approach the critical curvilinearity but not exceed it. The values of the coefficients are picked with reserve to provide stability for long runs.

Figure 2 Probability of the curvilinearity η_ξξ. The thick curve is calculated with a full 3-D model; the thin curve is the probability calculated over an ensemble of linear modes with the same spectrum.

We do not think that the suggested breaking parameterization is a final solution to the problem. Other schemes will be tested in the next version of the model. However, the results presented below show that the scheme is reliable and provides a realistic energy dissipation rate.

4 Calculations and results

The elevation and surface velocity potential fields are approximated in the current calculations by M_x=256 and M_y=128 modes in directions x and y. The corresponding grid includes $N_x\times N_y=(\mathrm{1024}\times \mathrm{512})$ knots. The vertical derivatives are approximated at a vertical stretched grid $d{\mathit{\zeta }}_{\mathrm{j}+\mathrm{1}}=\mathit{\chi }d{\mathit{\zeta }}_{\mathrm{j}},(j=\mathrm{1},\mathrm{2},\mathrm{3}\mathrm{\dots },L_w)$, where ν=1.2 and L_w=10. A small number of levels used for the solution of the equation for a nonlinear component of the velocity potential are possible because just a surface vertical derivative for the velocity potential $\partial \mathrm{\Phi }/\partial \mathit{\zeta }\left(\mathit{\zeta }=\mathrm{0}\right)$ is required. The velocity potential mainly consists of an analytical component $\overline{\mathit{\phi }}$, while a nonlinear component provides only a small correction. To reach the accuracy of the solution $\mathit{\epsilon }={\mathrm{10}}^{-\mathrm{6}}$ for Eq. (11), no more than two iterations were usually sufficient.

The parameters chosen were used for solution of a problem of a wave field evolution over the acceptable time (of the order of 10 days). The initial conditions were assigned on the basis of the empirical spectrum JONSWAP (Hasselmann et al., 1973) with a maximum placed at the wave number k_p=100 with the angle spreading (coshψ)²⁵⁶. The details of the initial conditions are of no importance because an initial energy level is quite low.

The total energy of a wave motion $E=E_{\mathrm{p}}+E_{\mathrm{k}}$ (E_p is a potential energy, while E_k is a kinetic energy) is calculated with the following formulas:

$$\begin{array}{}\text{(28)}& E_{\mathrm{p}}=\mathrm{0.25}\stackrel{\mathrm{‾}}{{\mathit{\eta }}^{\mathrm{2}}},E_{\mathrm{k}}=\mathrm{0.5}\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\left({\mathit{\phi }}_x^{\mathrm{2}}+{\mathit{\phi }}_y^{\mathrm{2}}+{\mathit{\phi }}_z^{\mathrm{2}}\right)}},\end{array}$$

where a single bar denotes averaging over the ξ and ϑ coordinates, and a double bar denotes averaging over the entire volume. The derivatives in Eq. (25) are calculated according to the transformation (1). An equation of the integral energy $E=E_{\mathrm{p}}+E_{\mathrm{k}}$ evolution can be represented in the following form:

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}E}{\mathrm{d}t}}=\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{I}}+\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}}}+\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{t}}}}+\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{N}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{I}}$ is the integral input of energy from wind (Eqs. 14–18); $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}}}$ is the rate of the energy dissipation due to the wave breaking (Eqs. 24–27); $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{t}}}}$ is the rate of the energy dissipation due to filtration of high-wave-number modes (tail dissipation, Eqs. 19–23); $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{N}}$ is an 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 form for calculation of the energy transformation 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 $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{I}}$ is calculated by the following relation:

$$\begin{array}{}\text{(30)}& \stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{I}}={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }t}}\left(\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{E^{t+\mathrm{\Delta }t}}}-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{E^t}}\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{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). For calculation of the dissipation rate due to filtration, the right-hand side of the equations contains just the terms introduced in Eqs. (19)–(23), while for calculation of the effects of breaking, only the terms introduced in Eqs. (24)–(27) are in use.

An evolution of the characteristics calculated by Eq. (30) is shown in Fig. 3. The sharp variations in all the characteristics at t<50 can probably be explained by adjustment of the initial linear fields to the nonlinearity. Up to the end of integration, the sum of all energy transition terms (tail dissipation $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{t}}}}$, breaking dissipation $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}}}$ and energy input $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{I}})$ approaches zero (curve 4), and the energy growth E (curve 5) stops. Then the energy tends to decrease, but we are not sure about the nature of this effect. Such behavior can be explained by a fluctuating character of mutual adjustment of input and dissipation or simply by deterioration of the approximation because of the downshifting process. Note that opposite to a more or less monotonic behavior of the tail dissipation (curve 1), the breaking dissipation is highly intermittent, which is consistent with the common views on the wave breaking nature.

The data on the evolution of the wave spectrum are shown in Fig. 4. A 2-D wave spectrum S(k,l) $\left(\mathrm{0}\le k\le M_x,-M_y\le l\le M_y\right)$ averaged over 13 time intervals of length equal to Δt≈100 was transferred to the polar coordinates $S_p(\mathit{\psi },r)\left(-\mathit{\pi }/\mathrm{2}\le \mathit{\psi }\le \mathit{\pi }/\mathrm{2},\mathrm{0}\le r\le M_x\right)$ and then averaged over the angle ψ to obtain a 1-D spectrum S_h(r):

$$\begin{array}{}\text{(31)}& S_h\left(r\right)=\sum S_p\left(\mathit{\psi },r\right)r\mathrm{\Delta }\mathit{\psi }.\end{array}$$

An angle ψ=0 coincides with the direction of wind U, $\mathrm{\Delta }\mathit{\psi }=\mathit{\pi }/\mathrm{180}.$

Figure 3 Evolution of integral characteristics of the solution, a rate of evolution of the integral energy multiplied by 10⁷ due to 1 – tail dissipation D_t (Eqs. 19–23), 2 – breaking dissipation D_b (Eqs. 24–27), 3 – input of energy from wind I (Eqs. 14–18) and 4 – balance of energy $I+D_{\mathrm{t}}+D_{\mathrm{b}}$. Curve 5 shows the evolution of wave energy 10⁵E. Grey vertical bars show the instantaneous values; the thick curve shows the smoothed behavior.

Figure 4 The wave spectra S_h(r) integrated over angle ψ in the polar coordinates and averaged over the consequent intervals of length for about 100 units of the nondimensional time t. The spectra grow and shift from right to left.

As seen, each spectrum consists of separated peaks and holes¹. This phenomenon was first observed and discussed by Chalikov et al. (2014). The repeated calculations with different resolution showed that such a structure of the 2-D spectrum is typical. It cannot be explained by a 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 at different locations in a Fourier space.

Another presentation is given in Fig. 6 in which the log₁₀(S(ψ,r)), averaged over the successive seven-period length Δt=200, is given. The first panel with a mark of 0 refers to the initial conditions. The disturbances within the range ($\mathrm{125}<k<\mathrm{150}$) reflect the initial adjustment of input and dissipation at a high-wave-number slope of spectrum. The pictures characterize the downshifting and angle spreading of the spectrum well due to the nonlinear interactions.

Figure 5 Sequence of 3-D images of lg₁₀(S(k,l)), in which each panel corresponds to a single curve in Fig. 3. The left side refers to the wave number $l\left(-M_y\le l\le M_y\right)$ and the front side refers to $k\left(\mathrm{0}\le k\le M\right)$. The numbers indicate the end of the time interval expressed in hundreds of nondimensional time units.

Figure 6 Sequence of 2-D images of lg₁₀(S(k,l)) averaged over the consequent seven periods of length Δt=200. The numbers indicate the period of averaging (first panel marked 0, refers to the initial conditions). The horizontal and vertical axes correspond to the wave numbers k and l, respectively.

Evolution of the integrated-over-angles ψ wave spectrum S_h(r) can be described with the equation

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\mathrm{d}{S}_h\left(r\right)}{\mathrm{d}t}}=I\left(r\right)+D_{\mathrm{t}}\left(r\right)+D_{\mathrm{b}}\left(r\right)+N\left(r\right),\end{array}$$

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

$$\begin{array}{}\text{(33)}& I\left(k,l\right)=\left(S_{\mathrm{c}}^{t+\mathrm{\Delta }t}\left(k,l\right)-S_{\mathrm{c}}^t\left(k,l\right)\right)/\mathrm{\Delta }t,\end{array}$$

where S_c(k_x,k_y) is the spectrum of the columnar energy calculated by the relation

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{S}_{\mathrm{c}}(k,l)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(h_{k,l}^{\mathrm{2}}+h_{-k,-l}^{\mathrm{2}}+\underset{-H}{\overset{\mathrm{0}}{\int }}.\\ \text{(34)}& {\displaystyle }& {\displaystyle }\left.(u_{k,l}^{\mathrm{2}}+u_{-k,-l}^{\mathrm{2}}+v_{k,l}^{\mathrm{2}}+v_{-k,-l}^{\mathrm{2}}+w_{k,l}^{\mathrm{2}}+w_{-k,-l}^{\mathrm{2}})d\mathit{\zeta }\right),\end{array}$$

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

$$\begin{array}{}\text{(35)}& u={\mathit{\phi }}_{\mathit{\xi }}+{\mathit{\phi }}_{\mathit{\zeta }}{\mathit{\eta }}_{\mathit{\xi }},v={\mathit{\phi }}_{\mathit{\vartheta }}+{\mathit{\phi }}_{\mathit{\zeta }}{\mathit{\eta }}_{\mathit{\vartheta }},w={\mathit{\phi }}_{\mathit{\zeta }},\end{array}$$

and $u_{k,l},v_{k,l}\mathrm{and}{w}_{k,l}$ are their Fourier coefficients.

For calculation of I(k,l) the fictitious time steps Δt are made only with a term responsible for the energy input, i.e., surface pressure p. A spectrum I(k,l) was averaged over the periods Δt≈100, then transformed into a polar coordinate system and integrated in a Fourier space over angles ψ within the interval $\left(-\mathit{\pi }/\mathrm{2},\mathit{\pi }/\mathrm{2}\right).$

Evolution of the input spectra (Fig. 7) is in general similar to that of the wave spectra shown in Fig. 4. Note that the maximum of the spectra is located at the maximum of the wave spectra since the input depends mainly on the spectral density, while the dependence on frequency is less important.

An algorithm (Eq. 30) was applied for calculation of the dissipation spectra due to dumping of a high-wave-number part of spectrum (tail dissipation) and for calculation of the spectrum of the breaking dissipation. In the first case, the fictitious time step was made taking into account the terms described by Eqs. (19)–(23), while in the second case the time step was made using the terms described by Eqs. (24)–(27).

The spectra of the tail dissipation calculated similarly to the spectra I(r) are shown in Fig. 8. The dissipation occurs at the periphery of the spectrum, outside an ellipse with semiaxes d_mM_x and d_mM_y². This is why such dissipation, averaged over angles, seems to affect the middle part of a 1-D spectrum. The tail dissipation effectively stabilizes the solution.

Figure 7 The spectrum of energy input I(r) integrated over angle ψ in the polar coordinates and averaged over the consequent intervals of length for about 100 units of the nondimensional time t.

Figure 8 The tail dissipation spectra D_t(r) integrated over angle ψ in the polar coordinates and averaged over the consequent intervals of length for about 100 units of the nondimensional time t.

The breaking dissipation averaged over angles is presented in Fig. 8. As seen, the breaking dissipation has a maximum at the spectral peak. This does not mean that in the vicinity of the wave peak the probability of large curvilinearity is quite high. A high rate of the breaking dissipation can be explained by high wave energy in the vicinity of the wave peak. The energy lost through the breaking, described by the diffusion mechanism, correlates with the energy of breaking waves. Opposite to the high-wave-number dissipation which regulates the shape of the spectral tail, the breaking dissipation forms the main energy-containing part of the spectrum.

The diffusion mechanism suggested in Eqs. (24)–(27) modifies an elevation and surface stream function in close vicinity of the breaking point. The amplitudes of side perturbation are small and decrease very quickly throughout the distance from the breaking point.

An example of the profile of the energy input due to the breaking D_b(x) is given in Fig. 10. As seen, the energy input fluctuates around the breaking point. A diffusion operator chosen for the breaking parameterization not only decreases total energy but also redistributes the energy between Fourier modes in a Fourier space.

Figure 9 The breaking dissipation spectra D_b(r) integrated over angle ψ in the polar coordinates and averaged over the consequent intervals of length for about 100 units of the nondimensional time t.

In general, for the specific conditions considered in this paper, the breaking is an occasional process taking place in a small part of the domain. The kurtosis of the input energy due to the breaking D_b(ξ,ϑ), i.e., the value

$$\begin{array}{}\text{(36)}& \mathrm{Ku}=\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}^{\mathrm{4}}}}{\left(\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{D_{\mathrm{b}}^{\mathrm{2}}}}\right)}^{-\mathrm{2}}-\mathrm{3},\end{array}$$

is of the order of 10³, which corresponds to a plain function with occasional separated peaks.

The number of breaking points in terms of percentage of the total number of points is given in Fig. 11. As seen, the number of breaking events decreases to t=600 and then increases till the end of the calculations. The number of breaking events is not directly connected with the intensity of breaking, which is seen when comparing Fig. 11 and curve 2 in Fig. 3.

An integral term describing a nonlinear interaction $\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{N}}$ in Eq. (29) is small (compared with the local values of N_k,l), but the magnitude of spectrum N(r) is comparable with input I(r) and dissipation terms D_t(r) and D_b(r). The presentation of term N(r) in the form shown in Figs. 6–8 is not clear. This is why the spectra 10⁸N(r) averaged over interval Δt=100 are plotted separately in Fig. 11 for the last eight intervals (thick curves) together with a wave spectrum 10⁶S_h(r). In general, the shapes of spectrum N(r) agree with the conclusions of the quasi-linear Hasselmann (1962) theory (Hasselmann et al., 1985). At a low-wave-number slope of spectrum the nonlinear influx of energy is positive, while at the opposite slope it is negative. This process produces the shifting of the spectrum to a lower wave number (downshifting). Opposite to the Hasselmann's theory, these results are obtained by solution of the full three-dimensional equations. It would be interesting to compare our results with the calculations of Hasselmann's integral. Unfortunately, neither of the existing programs of such a type permits calculations with the high resolution that was used in the current model. Note that the nonlinear interactions also produce widening of the spectrum.

Figure 10 Example of the energy input due to the breaking D_b(x).

As can be seen, the nonlinearity is quite an important property of surface waves. The contribution of nonlinearity can be estimated, for example, by comparison of the kinetic energy of a linear component $E_l=\mathrm{0.5}\stackrel{\mathrm{‾}}{\stackrel{\mathrm{‾}}{\left({\overline{\mathit{\phi }}}_x^{\mathrm{2}}+{\overline{\mathit{\phi }}}_y^{\mathrm{2}}+{\overline{\mathit{\phi }}}_z^{\mathrm{2}}\right)}}$ and the total kinetic energy E_k (Fig. 13). A ratio E_l∕E_k as a function of time remains very close to 1, which proves that the nonlinear part of energy makes up just a small percentage of the total energy. It does not mean that the role of the nonlinearity is small; its influence can manifest itself over large timescales.

The time evolution of the integral spectral characteristics is presented in Fig. 14. Curve 1 corresponds to the weighted frequency ω_w

$$\begin{array}{}\text{(37)}& {\mathit{\omega }}_w={\left({\displaystyle \frac{\int kS\mathrm{d}k\mathrm{d}l}{\int S\mathrm{d}k\mathrm{d}l}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where integrals are taken over the entire Fourier domain. The value ω_w is not sensitive to the details of the spectrum; hence, it characterizes the position of spectrum and its shifting well. Curve 2 describes an evolution of the spectral maximum. The step shape of the curve corresponds to the fundamental property of downshifting. Opposite to the common views, the development of spectrum occurs not monotonically but by the 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 a spectral model (Rogers et al., 2012). Curve 3 describes the change of total energy $E=E_{\mathrm{p}}+E_{\mathrm{k}}$. As seen, all three curves have a tendency to slow down the evolution rate. Then the energy tends to decrease, but we are not sure about the nature of this effect. Such behavior can be explained by a fluctuating character of mutual adjustment of input and dissipation or simply by deterioration of the approximation because of the downshifting process. 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 a JONSWAP experiment. Despite large scatter, the data allow us to construct empirical approximations of a wave spectrum, as well as to investigate the evolution of a spectrum as a function of fetch F. In particular, it is suggested that the frequency of a spectral peak changes as $F^{-\mathrm{1}/\mathrm{3}}$, while the full energy grows linearly with F. Neither of the dependences can be exact since they do not take into account the approach to a stationary regime. In addition, the dependence of frequency on fetch is singular at F=0.

Figure 11 Evolution of a number of the wave breaking events N_b expressed in percentage of the number of grid points N_x×N_y.

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

$$\begin{array}{}\text{(38)}& F=\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{c}_{\mathrm{p}}\mathrm{d}t\end{array}$$

The JONSWAP dependencies for the frequency of a spectral peak ω_p and the full energy E are shown in Fig. 14 by thin curves. Dependence ${\mathit{\omega }}_{\mathrm{p}}\sim F^{\mathrm{1}/\mathrm{3}}$ is qualitatively valid. A dependence of the total energy on fetch does not look like a linear one, but it is worth noting that the JONSWAP dependence is evidently inapplicable to a very small and large fetch.

Figure 12 Sequence of wave spectra S_h(r) (thick curves) and a nonlinear input term N(r) (thin curves) averaged over the eight consequent intervals of length Δt=100 starting from the sixth period.

Figure 13 Time evolution of the ratio E_l∕E_k.

Figure 14 Time evolution of weighted frequency ω_w (1) (Eq. 34), the spectral peak frequency ${\mathit{\omega }}_{\mathrm{p}}=k_{\mathrm{p}}^{\mathrm{1}/\mathrm{2}}$ (2) and full energy E(3) (Eq. 28). Thin curves correspond to the empirical dependences for the peak wave number and energy. F is a distance passed by the spectral peak.

5 Discussion

A model based on the full three-dimensional equations of a potential motion with a free surface was used for simulation of development of wave fields. The model is written in a surface-following nonstationary nonorthogonal coordinate system. The details of a numerical scheme and the results of the model validation were described in Chalikov et al. (2014). The main difference between the given model and HOS model (Ducroset et al., 2016) is that our model is based on a direct solution of the 3-D equations for the velocity potential. This approach is similar to that developed at the Technical University of Denmark (TUD; see Engsig-Karup et al., 2009). Actually, the models developed at TUD are targeted at the solution of a variety of problems including such problems as the modeling of wave interaction with submerged objects and the simulation of a wave regime in basins with real shape and topography.

In the current paper a three-dimensional model was used for simulation of the development of a wave field under the action of wind and dissipation. The input energy is described by a single term, i.e., surface pressure p in Eq. (4). It is traditionally assumed that the complex pressure amplitude in a Fourier space is linearly connected with the complex elevation amplitude through a complex coefficient called the β function. Such simple formulations can be imperfect. First, it is assumed that the wave field is represented by a superposition of linear modes with the slowly changing amplitudes and the phase velocity obeying a linear dispersive relation. This assumption is valid only for a low-frequency part of the spectrum. In reality, the amplitudes of the medium- and high-frequency modes undergo fluctuations created by reversible interactions. A solid dispersion relation does not connect their phase velocities with a wave number. In addition, it is also quite possible that a suggestion of the linearity of the connection between the pressure and elevation amplitudes is not precise, i.e., the β function can depend on the amplitudes of modes.

We are not familiar with any observational data that can be used for the formulation of a statistically provided scheme for calculation of the input energy to waves. The only method that can give more or less reliable results is the mathematical modeling of the statistical structure of a turbulent boundary layer above a curvilinear moving surface whose characteristics satisfy the kinematic conditions. The method described above is based on several millions of values of the pressure referred strictly to the surface. As a whole, the problem of a boundary layer seems even more complicated than the wave problem itself. Some early attempts to solve this problem were made on the basis of a finite-difference two-dimensional model of a boundary layer written in a simple surface following the coordinate (see review Chalikov, 1986). Waves were assigned as a superposition of linear modes with random phases, corresponding to the empirical wave spectrum. This approach was not quite accurate since it did not take into account the nonlinear properties of surface (for example, the sharpness of real waves and the absence of a dispersive relation for the waves of medium and high frequencies). The next step was the formulation of coupled models for a boundary layer and potential waves, both written in the conformal coordinates (Chalikov and Rainchik, 2010). The calculations showed that the pressure field consists mostly of random fluctuations not directly connected with the waves. A small part of these fluctuations are in phase with the surface disturbances. The calculated values of β in Eq. (13) have large dispersion. However, since the volume of data was very large, the shape of the β function was found with a high level of accuracy. Probably, the approximation of β used in the current work can be considered most adequate. We are planning additional investigations based on the coupled wind–wave models. The next step in the investigations of wave boundary layer (WBL) should use a three-dimensional LES approach. Note that even the availability of a large volume of data on the structure of WBL does not make the problem of parameterization of wind input in the spectral wave models easily solvable since the pressure is characterized by a broad continuous spectrum created by nonlinearity.

The wave breaking is obviously even more complicated than the input energy. Nevertheless, this problem can be simplified, if the common ideas used in the numerical fluid mechanics are accepted. For example, in the LES modeling a more or less artificial viscosity is introduced to prevent too large local velocity gradients. In fact the numerical instability terminating computations precedes the wave breaking. Hence, the scheme should prevent the breaking approach to preserve stability of the numerical scheme. Hence, a wave model should contain the algorithms preventing the appearance of too large slopes. A criterion of breaking is introduced not for recognizing the breaking itself, but for the choice of places where it might happen (or, unfortunately, might not happen). Finally, the algorithm should produce the local smoothing of elevation (and the surface potential). The algorithm should be highly selective so that the “breaking” could occur within narrow intervals and not affect the entire area. The exact criteria of the breaking events (most evident of them is the breaking itself) cannot be used for parameterization of breaking since in a coordinate system (1) the numerical instability occurs long before the breaking. In our opinion, the most sensitive parameter indicating potential instability is the curvilinearity (second derivative) of elevation.

In the current work, the breaking is parameterized by a diffusion algorithm with a nonlinear coefficient of diffusion providing high selectivity of the smoothing. We admit that such an approach can be realized in many different forms. The same situation is observed in a problem of the turbulence modeling for parameterization of subgrid scales. Note that the breaking dissipation in phase-resolving models is included in a more realistic manner than in spectral models. For example, the breaking is simulated in a physical space, which allows us to reduce the height and energy of the nonlinear waves composed of several modes. In spectral models the dissipation is distributed more or less arbitrarily over the entire spectrum. The spectral models sometimes include additional dissipation of short waves due to their modulation by long waves (Young and Babanin, 2006; Babanin et al., 2010). In the phase-resolving models this process has been included explicitly.

We can finally conclude that the physics included in wave models still rests on shaky ground. Nevertheless, the result of the calculations looks quite realistic, which convinces us that the approach deserves further development.

The numerical models of waves similar to that considered in this paper have a lot of important applications. First, they are a perfect tool for the development of physical parameterization schemes in spectral wave models. Second, a direct model can be used in future for the numerical simulation of wave processes in the basins of small and medium size. These investigations can be based on the HOS model (Ducrozet et al., 2016) or the model used in the current paper. However, the most universal approach seems to be developed at the Technical University of Denmark (see Engsig-Karup, 2009). Any model used for the long-term simulation of wave field evolution should include the algorithms describing transformation of energy similar to those considered in the current paper.