**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.

*High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind DOI: http://dx.doi.org/10.5772/intechopen.92262*

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},\tag{21}$$

$$\frac{\partial \rho\_{k,l}}{\partial \tau} = F\_{k,l} - \mu\_{k,l} \rho\_{k,l} \tag{22}$$

(where *Ek*,*<sup>l</sup>* and *Fk*,*<sup>l</sup>* are the Fourier amplitudes of the right-hand sides of Eqs. (4) and (5); the value of *μ<sup>k</sup>*,*<sup>l</sup>* is equal to zero inside the ellipse with semiaxes *dmMx* and *dmМy*; then, it grows quadratically with j j *k* up to the value *cm* and is equal to *cm* 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–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 \theta} B\_{\theta} \frac{\partial \eta}{\partial \theta} \right), \tag{23}$$

$$
\rho\_{\tau} = F\_{\rho} + f^{-1} \left( \frac{\partial}{\partial \xi} B\_{\xi} \frac{\partial \rho}{\partial \xi} + \frac{\partial}{\partial \theta} B\_{\theta} \frac{\partial \rho}{\partial \theta} \right), \tag{24}
$$

where *F<sup>η</sup>* and *F<sup>φ</sup>* are the right-hand sides of Eqs. (4) and (5) including the tail dissipation terms; *B<sup>ξ</sup>* and *B<sup>ϑ</sup>* 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<sup>ξ</sup>* and *B<sup>ϑ</sup>* depend nonlinearly on the curvilinearity

$$B\_{\xi} = \begin{cases} \mathcal{C}\_{B} \eta^{2}\_{\xi\xi} & \eta\_{\xi\xi} < \eta^{cr}\_{\xi\xi} \\ & \mathbf{0} & \eta\_{\xi\xi} \ge \eta^{cr}\_{\xi\xi} \end{cases} \tag{25}$$

$$B\_{\theta} = \begin{cases} \mathcal{C}\_{B} \eta^{2}\_{\theta\theta} & \eta\_{\theta\theta} < \eta^{cr}\_{\xi\xi} \\ & \mathbf{0} & \eta\_{\theta\theta} \ge \eta^{cr}\_{\xi\xi} \end{cases} \tag{26}$$

where the coefficients at *CB* <sup>¼</sup> <sup>0</sup>*:*05, *<sup>η</sup>cr ξξ* <sup>¼</sup> *<sup>η</sup>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.
