2. Derivation of a Simplified 2D Wave Model (HWM) for 2D Waves
The equations are written in the non-stationary surface-following non-
orthogonal coordinate system:
where
is a moving periodic wave surface defined by Fourier series
where
k and
l are the components of the vector wave number
,
are Fourier amplitudes for elevations
,
and
are the numbers of modes in the directions
and
, respectively, whereas
are the Fourier basic functions represented as a matrix:
The surface conditions for potential waves in a system of coordinates (1) at
take the following form:
where
is the velocity potential. Laplace equation for
at
turns into the Poisson equation
where
is the operator:
Strictly speaking, Equations (4)–(6) should be solved jointly, but since the Fourier amplitudes for elevation and for the surface velocity potential change slowly, some sort of time-splitting scheme is used: Equations (4) and (5) are used as evolutionary equations for and , whereas Equations (6) and (7) are used to satisfy the continuity conditions. It is accepted that in the process of the solution of Equation (7) the coefficients depending on are constants. Introducing the correction of and s with Equations (4) and (5) in the course of solution does not provide any advantages.
Equations (4)–(7) are written in a non-dimensional form by using the following scales: length L, where 2πL is a dimensional period in a horizontal direction; time (g is the acceleration of gravity) and the velocity potential . The pressure is normalized by water density so that the pressure scale is Lg. Equations (4)–(6) are self-similar to transformation with respect to L. The dimensional size of the domain is . All of the results presented in this paper are nondimensional.
It is suggested in [
7] that it is convenient to represent the velocity potential as a sum of two components: an analytical (‘linear’) component
and an arbitrary (‘non-linear’) component
:
The analytical component
satisfies Laplace equation:
with the known solution:
where
,
are Fourier coefficients for the surface potential
at
. The 3D solution of (9) satisfies the following boundary conditions for the surface and at infinite depth:
The presentation (8) is not used for solution of evolutionary Equations (4) and (5), because it does not provide any noticeable improvements of accuracy and speed of the calculations.
The nonlinear component satisfies an equation:
Note that the full potential
on the surface is used as a boundary condition for a linear component, so, a nonlinear component
on the surface is equal to 0. Thus, Equation (12) is solved with the boundary conditions:
Hence, on the surface
Equation (12) takes the form:
The derivatives of a linear component
in (7) are calculated analytically. A method of solution combines a 2D Fourier transform method in the ‘horizontal surfaces’ and a second-order finite-difference approximation on the stretched staggered vertical grid providing high accuracy in the vicinity of the surface. An Equation (12) in FWM is solved as Poisson equations with iterations over the right-hand side by the TDMA method [
10] with the prescribed accuracy. The number of iterations depends on a maximum of the local surface steepness, and for high steepness it does not exceed three.
A 3D model described above is potentially exact. Its only drawback is a low speed of performance. Below is a scheme that at a cost of minor loss of accuracy provides a significantly higher speed of calculation.
Since
, an Equation (12) for the velocity potential at
takes the form:
Here, the symbols
and
are used. It is remarkable that an Equation (15) is exact. From Equation (14) it follows that
hence, Equation (15) can be represented as follows:
where the term
depends only on a linear component
:
which is calculated analytically using Fourier presentation (10) by the formulas:
The presence of
in Equation (17) makes the system of equations unclosed. Previously, the author regarded a surface condition (17) as a means of validation of the numerical scheme for the 3D Full Wave Model (FWM) based on a numerical solution of an equation for the velocity potential (see example of such calculation in [
8]). This validation is not trivial since
is not a direct product of the Poisson equation solution. The finite-difference derivatives
and
after solution of the 3D equation for a nonlinear component of the velocity potential were calculated as follows:
Here, , is depth where is located, , the thickness of layer grows exponentially with depth : , whereas the value of is calculated from the condition where is the number of levels, and is depth of the domain ( is a wave number of a mode with the largest amplitude).
Since a relation (17) formally follows from the Poisson equation, the agreement between these values characterizes accuracy of the simplified scheme. The agreement cannot be very precise since the variables
and
are not a direct product of the solution, and they are calculated by the approximations (23) and (24) following the solution of the Poisson equation. Moreover, the approximations for solution of the Poisson equation and calculation (23) and (24) are different. However, the calculations of
on the basis of Poisson Equation (12) with the approximations (23) and (24) and calculations by a surface condition (17) showed a very good agreement (see Figure 1 in [
9]), which prompted us to carry out further analysis of Equation (18).
It was found that variables
and
turned out to be well connected to each other. Firstly, the simplest linear hypothesis was tested:
where
and
are constants. This dependence was considered in [
11] with the use of the 3D FWM model (described in [
1,
7,
8,
9]) that calculates the 3D structure of the velocity potential as well as 2D fields of
and
. The dependence between
and
can be studied in Fourier space as well as in a physical space. The results obtained in Fourier space turned out to be far less certain than in a physical space, i.e., the dependence between the Fourier coefficients for
and
showed considerable scatter. In our opinion, it happens because the medium and high wave numbers of Fourier amplitudes are generally unstable, whereas the verticals structure of a physical layer involved in orbital movement is universal because it is composed of a large number of modes.
The calculations performed in [
11] were repeated here at a higher grid resolution
and a number of vertical levels equal to 50. The dispersion of surface elevation
, as well as the integral steepness
s, variate in the ranges
and
, respectively. The fields included
pairs of values
and
. The statistical connection between those variables is shown in
Figure 1. It is important to underline that dependence (25) is valid for the entire wave field, i.e., the parameters
and
are close to constants. Otherwise, the whole approach based on a hypothesis (25) might be incorrect.
The data used for analysis included the 3D velocity potential fields and the values of and calculated by (23); the values of calculated by a relation (17) and some other characteristics (see below). Since a high value of for generation of vertical grid was used, the accuracy of calculation of the derivatives (23) and (24) was very high. The comparison with the analytical values of the derivatives shows that the errors of approximation (23) and (24) do not exceed .
The surface values of and are defined by an asymptotic behavior of a nonlinear component of the velocity potential at . Hence, the connection between them can be found in the close vicinity of .
A dimensional form of (25) is
where
and
are dimensional analogues of
and
, whereas
is the external scale of the problem. Formally, for any specific wave field which is more or less horizontally homogeneous, the dependence (25) is correct at appropriate choice of coefficients. However, the problem is that for different wave fields (for example, for different stages of wave field development) the variability of those coefficients can be significant. It happens because the scaling with external scale
does not reflect the internal properties of wave field. To provide a more universal approximation, it is necessary to introduce the scaling in terms of integral characteristics of a given wave field with no use of
.
For evaluation of connection between
and
, ten thousand short-term numerical experiments with FWM were performed. The Fourier resolution is equal to
modes; the number of levels for solution of Poisson equation equals 30, whereas a stretching coefficient
. The initial conditions for elevation were calculated on the basis of JONSWAP spectrum [
12] for the inverse wave age
(
is wind velocity,
is the phase velocity in peak of spectrum). The calculations were carried out for different values of peak wave numbers in the range
. The data collected are characterized by the nondimensional dispersion of elevation
in the range
and by steepness
in the range
and by dispersion of Laplacian
in the range
. The calculations performed with such wide variations of integral parameters showed that the values
variate in the range
, whereas the dependence of
on
is not linear as assumed in [
11]. Hence, the problem is reduced to finding the dependence
In fact, several additional parameters were tested, whereas (Equation (27)) and (Equation (29)) only turned out to be well connected with the ratio .
For generalization of the linear dependence (25) it is necessary to choose at least two parameters with a dimension of length taking into account internal characteristics of a wave field type (27)–(29). Finally, the best results were obtained in the following form
where
is a parameter
and function
is approximated by a formula
where
The function
is shown in
Figure 2.
An approximation (33) is also correct for dimensional variables, because it is independent of the external scale . The form of relation (31) and (32) and the constants in (33) were found on the basis of numerous numerical experiments with FWM (Equations (4)–(6)) and empirical selection of nondimensional variables, as well as functions and numerical parameters. That is why such simplified model is called the Heuristic Wave model (HWM).
The two-dimensional model, which is presumably able to replace the full 3D model, (4), (5), and (12), looks as follows:
Note that the right-hand side of Equation (36) contains full vertical velocity
as well as a linear component of vertical velocity
. An equation is represented in a form that is convenient for iterations. The iterations were performed until the next condition was reached.
where
is the right-hand side of Equation (36), whereas
is the number of iterations. In the majority of cases, the number of iterations was equal to two and never exceeded four. Note that fulfilment of condition (37) can be regulated by an appropriate choice of parameters in the breaking algorithm (see below).
A comparison of function
calculated in the course of simulation by FWM with the results of calculations by Formula (33) is given in
Figure 3. The agreement between these functions is quite satisfactory: a correlation coefficient between the model values
and those calculated by Equations (31)–(33) is equal to 0.989, the root-mean-square difference between them equal to 0.014.
In fact, the research described in this paper was initiated due to the fact discovered in the first version of the model [
11]. It was found that a nonlinear component of vertical velocity
(which is normally taken from solution of the Poisson equation) is a small value as compared with a full component
. However, the evolutionary Equations (34) and (35) contain a full component
only. It is illustrated in
Figure 4 by comparison of the cumulative probability both for the full component
and the nonlinear components
(in
Figure 4 they are marked as
and
, respectively). On the average,
is by one decimal order smaller than
. The main idea of the current approach is to replace the cumbersome algorithm for calculation of the vertical velocity equation
with a simple 2D Equation (36). Such trick reduces the total 3D problem to the 2D one, which saves us from the endless problems associated with solving Poisson equations numerically, as well as from any other ways of reproduction of 3D structure of the velocity potential. It allows us to create a model running much faster than the 3D model.
The most convincing is the comparison of two independent methods of calculation of vertical velocity
, i.e., a method based on solution of Poisson Equation (12) (marked as
) and another one—based on the surface condition (36) and Equations (31)–(33) (marked as
). That comparison is given in
Figure 5. The number of pairs of
and
used is
. The root-mean-square difference between
and
is equal to 0.0017.
The question is how well these simplified equations reproduce wave dynamics.
A simplified model is developed on the basis of the ‘exact’ model, both models having a similar structure and a set of parameters. The evolutionary equations of both models (4), (5), (34) and (35) are essentially the same. The difference between the models is determined by the simplified calculation of small nonlinear correction to the full vertical velocity introduced by Equations (31) and (33). A straightforward way of validation of the simplified model is performing the runs with the identical setting and initial conditions. A simple case of a quasi-stationary regime was considered.
The similarity of the models can be illustrated by difference
(
t is time)
and a correlation coefficient
between the fields of elevation generated by a full model
where
and
are dispersions of the fields
and
; the over-line denotes the averaging over the entire field. The right hand-sides of (38), (39) are calculated with the time interval equal to
. The evolution of
and
given in
Figure 6 shows that the fields preserve their closeness well up to
(10,000 time steps). Such a result turned out to be quite unexpected because such a close agreement was observed for large wave fields composed of thousands of nonstationary and dispersing Fourier modes. The comparison of FWM with a simplified model was also made in [
11],
Figure 5. Both models used the spectral resolution of
Fourier modes and
knots. In the initial conditions JONSWAP spectrum with the peak wave number
is assigned. It was obtained that the evolution of a correlation coefficient and a root-mean-square difference between the elevation fields was quite similar to those shown in
Figure 6.
Figure 6 proves that both models launched with the same boundary conditions produce a similar development for hundreds of wave peak periods. The advantage of a new scheme (Equations (31)–(33)) is its applicability to a wide range of external parameters.
The same comparison was made for FWM and HWM with a zero nonlinear correction to the vertical velocity , i.e., using no Equation (36). The run was quite stable but an agreement between FWM and such fastest version of a simplified model was observed at a two- or three-times shorter period. Essentially, it is not clear what is more responsible for the nonlinearity: either the presence of nonlinear terms in the evolutionary Equations (34) and (35), or introducing a nonlinear correction to the vertical velocity through the exact Equation (12) or a simplified algorithm ((31)–(33) and (36)).
3. An Example of Simulation of the Wave Field Development on the Basis of HFM
The comparison of the long-term simulations performed with FWM and HWM was made in a previous publication [
11]. That comparison proves that the accelerated and full models give the results quite similar to each other. The statistical characteristics of the solution were particularly well reproduced.
In the current work the calculations were performed with HWM model with the initiated algorithms for input and dissipation of energy. The algorithms for calculation of those terms are described in detail in [
1,
7,
8,
9]. The last minor modification of an algorithm for breaking is carried out in [
10]. Those algorithms are not described in the current paper.
The model uses
Fourier modes or
knots. The number of degrees of freedom (a minimum number of Fourier coefficients for elevation
and the surface velocity potential
) is 264,196. In the initial conditions the JONSWAP spectrum [
12] with a peak wave number
and the angle distribution proportional to
was assigned. It corresponds to nearly unidirected waves with very small energy. Note that the details of the initial conditions assigned at sufficiently high wave numbers are not important, since the model quickly develops its own wave spectrum.
The model, as well as FWM, is supplied with a developed system of algorithms for monitoring the solution and recording many functions in a spectral and grid form. The most important were the spectra of elevation and different terms describing input and dissipation of energy. An exact expression for total kinetic energy of wave motion
in the surface-following coordinates (1) can be represented in terms of surface variables
and
:
where a single bar denotes the averaging over the
and
coordinates. A relation (40) is easily derived with a use of Poisson Equation (7) and an expression for kinetic energy written in the coordinates (1). The total potential energy
is simply calculated by the formula
An equation of the integral energy
evolution can be represented in the following form:
where
is the integral input of energy from wind;
is a rate of energy dissipation due to wave breaking;
is a rate of energy dissipation due to filtration of the high-wave number modes (‘tail dissipation’);
is the integral effect of the nonlinear interactions described by the right-hand side of the equations when the surface pressure
is equal to zero. The differential forms for calculation of energy transformations can be, in principle, derived from Equations (4)–(6), but here a more convenient and simple method is applied. Different rates of the integral energy transformations can be calculated in the process of integration. For example, the value of
is calculated by the following relation:
where
is the integral energy of a wave field obtained after one time step with the right hand-side of Equation (35) containing only the surface pressure. This method is also applicable to the calculation of 1D spectra (i.e., the averaged over lateral wave number or along the direction in the polar coordinates) of different characteristics (see [
1]).
An evolution of integral characteristics as a function of time
(expressed in the number of peak wave period
in the initial conditions) is given in
Figure 7. The fastest change of all the characteristics occurs in the first
until the moment when the breaking dissipation
(curve 2) starts to develop. The input energy
(curve 1) remains nearly constant starting from
. The total balance
(curve 4) gradually decreases but it does not turn into zero, and the total energy
(curve 4) continues to grow. Curve (5) describes an evolution of the spectrum-weighted wave number
:
It is a more convenient characteristic of the position of the spectrum because it changes in time smoothly. Over the time of integration, decreases from the value to .
The nonlinear transformation of spectrum described by the right hand-sides of Equations (34) and (35) always produces the output of energy outside the computational domain in Fourier space. This process can be treated as dissipation. In our opinion, to keep the whole process under control, it is reasonable to maintain the total energy in the adiabatic part of the problem, governed by the initial Equations (4)–(6) with no input energy and dissipation. It can be performed by correction of all Fourier coefficients for the potential
and elevation
:
where a coefficient
is calculated by the formula
where
is full energy calculated by a set of the coefficients
before a time step and by
—after a time step. The correction (45) and (46) is performed at each time step before calculating and adding the input and dissipation terms. A typical value of the coefficient
is
. It is important that the correction (45) and (46) eliminates not only the loss of energy due to energy outflow from the computational domain but also due to the errors of approximations. A typical value of integral dissipation of this type was by around 2 decimal orders smaller than other energy transformation effects. Finally, the integral effect of nonlinear interactions was equal to zero (a thin straight line in
Figure 8).
Additional integral characteristics are shown in
Figure 8 where the computed values as a function of fetch
are compared with the same values calculated with a use of JONSWAP approximation [
12]. Since the periodic variables are used, there is no other alternative but to assume that the fetch is a distance passed by a peak of the wave spectrum:
The dependence of
on
can be obtained with a use of JONSWAP approximations:
The dependences (47) calculated by HWM and Eqauation (48) are shown in
Figure 8 by thick and thin curves, respectively. As seen, these curves are quite close to each other. The dependence of peak frequency (defined as the frequency in a maximum of wave spectrum) is shown by curves 2. A dashed thick curve is calculated with the HWM model, and a thin curve corresponds to the dependence following from JONSWAP approximation:
Opposite to Equation (44), the frequency
is a point value. Due to the irregularity of the 2D wave spectrum, the peak frequency changes in time not monotonically but in jumps, and sometimes
even moves back, i.e., to higher wave numbers. This is why the dependence
looks distorted. In
Figure 8 this curve is smoothed. However, an agreement with the experimental dependence is quite obvious. The JONSWAP approximation of dependence of energy
on fetch F is linear:
This dependence is shown in
Figure 8 by a dotted line 3, whereas the dependence calculated by EWM is shown by thick dots. Those dependences agree with each other only at early stages of wave development. Evidently, an approximation (50) cannot be correct for a long fetch because it does not take into account the saturation of the spectrum. The more reasonable form is:
where
is a function of the inverse wave age
(
is wind velocity).
In general, the results presented in
Figure 7 and
Figure 8 prove that an evolution of integral characteristics is reproduced satisfactorily.
The HWM also reproduces well the probability of the nondimensional surface elevation
. The height of large waves predicted by HWM is somewhat bigger than that predicted by FWM [
9] due to a larger total steepness of wave field and because the integration time was much longer. It is interesting to note that normalization of surface elevation by the significant wave height is so efficient that the probability distribution for the nondimensional surface elevation
turns out to be nearly universal, i.e., it does not depend on wave energy and just slightly depends on total wave steepness. Curves in
Figure 9 are calculated over
points; hence, it is not surprising that they are very smooth. In reality, the statistics of surface elevation shows very large scatter. The probability of the trough-to-crest height was calculated by a method of moving windows (see [
13,
14]). Each grey curve in
Figure 10 is calculated over individual wave field containing 524.288 points. The total number of curves is 594. As seen, the probability for different wave fields variates in a broad range: many of them do not contain ‘freak’ waves with the height of
at all, whereas in many cases their height can exceed such a high value as
and even reach a value
. This is why prediction of extreme waves for practical purposes is senseless with no indication of confidence intervals or distribution of probability for each value of
.
Figure 11 demonstrates statistical characteristics of positive and negative elevations as a function of the trough-to-crest height. As seen, the crest height
grows with the increase in total height
faster than the trough depth, which can be explained by nonlinearity of the process.
Figure 11 shows that on the average the trough-to-crest height
corresponds to the crest height
and the trough depth
. The dependences between those characteristics can be approximated by the formulas:
where tilde denotes normalization by the significant wave height
. Note that the first and third coefficients in approximation (52) coincide, which proves high quality of the ensemble. A dashed curve corresponds to the same data obtained with FWM [
10]. The difference between the new and old data is growing for large values of
. The old data contain 41.677 points, whereas a new array contains 1.264.441 points.
In a linear wave field (i.e., superposition of linear modes with random phases), the probabilities of crest height and trough depths are identical, but the probability of the trough-to-crest height is close to that for the nonlinear waves in a wave field with the same energy. This fact reliably proved with extensive calculations using the accurate two-dimensional and three-dimensional models [
13,
14,
15], has not been explained yet. This fact seems to be directly related to the physics of extreme waves. The nonlinearity manifests itself in vertical asymmetry of waves, but not in the probability of full height.
The next three figures show the 2D spectra obtained by averaging over six consecutive periods of length
(about 1600 initial peak period). Each spectrum is calculated by averaging 100 spectra transformed into the polar coordinates
. The surface elevation spectra
are given in
Figure 12. The colors correspond to different values of
where
is a maximum of the last of the six spectra. A minimum value (black color) corresponds to the values less than
. It is seen that the spectrum grows and approaches the point
. Despite the severe averaging, the spectra look highly irregular. They contain multiple peaks and holes (well seen in panels 5 and 6). The discussion of this effect is given in [
7,
13].
The spectra of the energy input shown in
Figure 13 calculated similar to wave spectrum are much wider than the wave spectrum, because the input of energy decreases with growth of wave number slower than the wave spectrum. It happens because
-function (see [
16]) depends on wave number quadratically. For large fetches the areas of negative wind input in low wave numbers (i.e., a flux of energy from fast wave modes to wind) usually arise, but it cannot be shown in the given figures.
The spectra of total dissipation rate equal to the sum of breaking and tail dissipation, are shown in
Figure 14. A maximum of dissipation always falls on a maximum of wave spectrum. However, dissipation, as well as the input energy, occurs everywhere including the spectral tail. Here, as in all our previous calculations, we did not reproduce an adiabatic regime at the spectral tail, which excludes any possibility to draw analogy to Kolmogorov’s spectrum. The attenuation of energy in the spectral tail depends on a shape of entire spectrum, as well as on wind velocity, wave age and other factors.
The process of downshifting is well seen in
Figure 15 where the evolution of 1D wave spectra (black curves) and the spectra of the nonlinear transformation rate (red curves) in the vicinity of a peak averaged over three consequent periods, each of
length, are given. All the spectra are calculated over 198 records of the two-dimensional spectra transformed into the polar coordinates, summated over angles and represented as a function of radius
. The locations of a wave spectrum maximum are indicated by thin vertical lines. Note that the analysis of the spectra for a nonlinear transformation rate is not simple, because the dispersion of reversible exchanges of energy between the modes is much higher than the rate of residual irreversible interaction. The accuracy of calculations grows with increase in spectral resolution. The recently launched calculations with HWF using spectral resolution of
modes gave the averaged spectral distribution much smoother than that in
Figure 15. In general, the spectra in
Figure 15 agree with the known regularities: the energy from a right slope of wave spectrum is transferred to the left slope, i.e., to longer waves. It qualitatively confirms Hasselmann’s theory [
17] considering the four-wave interactions. Opposite to that theory, actual interactions occur between all the spectral modes within the scope of the Fourier transform method. Unfortunately, the quantitative comparison with the theory is impossible since neither of the codes for calculation of Hasselmann integral can operate with such a high number of modes as used in phase-resolving models.
Most of the energy coming from wind to waves dissipates because of wave breaking. A currently used parameterization method for wave breaking is based on a diffusion algorithm with a variable coefficient of diffusion applied to both the elevation and surface velocity potential. The breaking occurs when Laplacian of elevation
becomes smaller than the negative value
(in the current model
= −80). A diffusion coefficient
is defined as follows:
where
is a horizontal step, and
is a coefficient. The parameters
and
are selected with consideration of a rate of growth of total energy. It is remarkable that the integral rate of breaking dissipation is not too sensitive to the choice of parameters because the effects of frequency of breaking (regulated by
) compensate for the intensity of breaking (depending on
). The rate of breaking seems to depend mostly on a rate of the energy input: most of the energy should be removed to maintain the numerical stability in the model, i.e., to prevent appearance of too sharp waves. Obviously, in a real wave field similar effects occur: dissipation absorbs extra energy obtained from the wind. Though a suggested algorithm seems very simple, the scheme performs all that is required.
Note that some part of the released horizontal momentum and energy should be transferred to horizontal motion (current). This effect cannot be taken into account in a potential model. The energy released after breaking is partially dissipated and partially distributed in space and between the wave modes, thus additionally contributing to the nonlinear interactions. The effects of breaking are presented in
Figure 16 where the number of breaking cases (% of the total number of knots) and the integral rate of energy loss at each time step are presented. As seen, both characteristics grow in time with the tendency for stabilization by the end of a run.
At an early stage of the work on parameterization of breaking the attempts were made to construct a scheme based on the geometrical and mechanical consideration with direct calculations of the stability criterion, a detaching volume of water and the possible consequences. Finally, we came to the conclusion that such a scheme would never work properly or be sufficiently robust to be used at millions of grid point and time steps. It is more reliable to adhere to the standard methods adopted in the geophysical fluid mechanics.
The width of spectrum can be characterized by a function
(see [
18,
19]).
where integrals have taken over the domain
. The wave spectra as the functions of frequency
normalized by peak frequency
for the same three periods shown in
Figure 15 are given. As seen, the
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, as well as a well-pronounced minimum in the spectral peak and a relatively slow growth above the spectral peak. The decrease in
at high frequencies is probably caused by the high-frequency dumping. The angle distribution was investigated in [
18,
19]. The approximations of
from the different sources collected in [
20] show considerable scatter, but the general features are quite similar to those shown in
Figure 17.
An example of a simulated elevation field
normalized by the significant wave height
is given in
Figure 18. In the initial conditions, small unidirected waves at inverse wave age
were assigned. In
Figure 18 we can see a developed three-dimensional field with the energy by two decimal orders higher than that in the initial conditions, and the broad angle distribution. Black dots mark the points where breaking occurs at a given time step.
4. Conclusions
The three-dimensional modeling of surface waves based on full nonlinear equations is a powerful tool for the investigation of wave processes, development of the parameterization schemes for the phase-resolving and spectral models and direct simulation of wave regimes in small basins. This type of modeling is rapidly developing. However, all approaches have a common limitation, i.e., low computational efficiency. Working with such models, even with modest resolution, turns into an endless waiting of the results. This property of the models slows down their improvement, in particular, development of the parameterization schemes for physical processes, since such work requires multiple repetitions of the runs. Such schemes significantly depend on resolution of a model, which limits the possibility of the research with low-resolution models.
The current paper continues to develop a new approach to the phase-resolving modeling of two-dimensional periodic wave fields at infinite depth. The basic concept of the scheme follows from presentation of the velocity potential as a sum of the linear and nonlinear components suggested in [
7]. The solution for a linear component is known; hence, a nonlinear component should be calculated through Poisson equation with a zero-boundary condition on the surface. It was observed in [
10] that such approach offers a new way to simplify the calculation by considering of 2D Poisson equation on the surface. The equation that can be treated as additional exact surface condition contains both the first
and second
vertical derivatives of the potential. Thus, the system of equations remains unclosed. It was empirically discovered that these variables are closely connected to each other. The linear dependence between
and
was checked in [
10]. It was shown that the use of such hypotheses leads to formulation of a closed system of equations.
In this paper, the above approach was improved by introducing a more flexible dependence of a ratio
on the nondimensional integral parameters of the wave field. The comparison of full vertical velocity
calculated using a new method, with that calculated by FWM, shows good agreement (
Figure 5). The point-to-point comparison of runs with HWM and FWM starting from the identical initial conditions proves that the solutions agree over thousands of time steps (
Figure 6). A long high-resolution run with HWM revealed that a simplified model produces the statistical results similar to those obtained with FWM [
9]. The main advantage of the new approach is high computational efficiency. The ratio of time spent for calculation with identical FWM and HWM variated in the range of 50–100. This ratio depends on complexity of parameterization of physics and volume/frequency of the processed and recorded information. Some acceleration can be reached by replacing a Runge–Kutta time integration scheme with a high accuracy one-step scheme.
The model suggested is not mathematically precise since it is based on the empirical (in computational sense) connection between the local variables, i.e., vertical velocity and its vertical derivative. The fact that those variables are connected by an equation containing horizontal derivatives indicates that the local connection cannot be exact. To discuss this problem, it is necessary to answer a general question: what does the concept of precision mean?
It should be noted that surface waves are often investigated on the basis of different simplified equations such as a linearized Navier–Stokes equation, Korteweg de Vrie equation, a linear and nonlinear Schrödinger (NLS) equation and others. No one actually knows which real properties of sea waves are omitted in those models. In some simplified approaches (e.g., in the models based on a nonlinear Schrödinger equation (see [
21]), the instability of high waves is missing; hence, the amplitudes of simulated waves can be unrealistically large. In reality, an excessive growth of waves is prevented due to the energy focusing and wave-breaking instability. The models based on the full potential equations predict the probability distribution for large waves quite satisfactorily. The main obstacle to a widespread use of the models based on complete equations is their low computational efficiency. An approach suggested here eliminates this obstacle by excluding reproduction of a three-dimensional flow structure. However, nothing prevents us from using a combined model with the on-demand calculations of a three-dimensional field of the velocity potential. A simplified model gives almost the same statistical results as a full model. Essentially, the difference between the statistical results given by 3D and 2D models are similar to those obtained with a 3D model initiated at different sets of random phases. It is easy to see that the equations suggested above are completely similar to the full equations, with one exception, i.e., a small non-linear correction to the total vertical velocity is calculated not from a 3D Poisson Equation (12) but from a simple 2D Equation (17). It can be added that a simplified model has a much simpler structure than FWM and is easily programmed. It is known that in most cases the data on 3D velocity potential are not required, but, if necessary, they can be easily obtained periodically in a process of model run or by the post-processing using the restart records and the corresponding subroutines adopted from 3D FWM.
It is obvious that an approach developed here is not fully precise. For example, it cannot be applied for individual cases with a small number of modes, for example, for simulation of a steep Stokes wave, as it was demonstrated for FWM in [
7,
9]. The 2D model is intended for simulation of the statistically homogeneous multi-mode wave field.
To continue the discussion on the accuracy problem, it should be noted that the full 3D adiabatic model is also inaccurate, because it is suitable for relatively short time intervals only. Real waves receive energy from wind and dissipate. Currently, such processes are poorly studied because they are evidently more complicated than the wave movement itself. Therefore, neither of the models can be considered as fully adequate. Especially for the followers of high accuracy it should be added that the full potential equations are far from being accurate even for the waves with no input energy and breaking. Real waves develop in a flow with the velocity shift, and they exchange energy with vortex movements [
22]. Therefore, they are better described by the Navier–Stokes equations. As a result, we obtain a variety of problems associated with the turbulence. It may turn out that the simplifications used in the construction of this model may be quite insignificant, as compared with other disadvantages.
Certainly, there are many ways to improve the model, first of all, for example, testing of alternative schemes for calculation of a nonlinear velocity component (Equations (31)–(33)). The attempts were made to construct the closure scheme for a surface Poisson equation in Fourier space. This way turned out to be unrealistic, since the local structure of the velocity potential depends not only on Fourier variables but also on a local full elevation.
The most evident advantage of a new model is the absence of calculation of the 3D structure of the velocity potential, which increases the speed of calculations by around two orders. This easily programmable model can be used to improve the physical components of the model and the long-term simulations of a multimode wave field. The model can be also used for acceleration of the calculations as a block of an exact 3D model (such as FWM). This coupling provides the possibility of runtime tuning of the parameters of the closure scheme for a surface Poisson equation.
The scope of practical applications of such model is quite wide. First of all, it should be used for building up of parametrization of physical processes in spectral and phase-resolving models. HWM, due to its high performance, can be used as a component of wave prediction models (such as WAM or WAVEWATCH) for interpretation of spectral information in terms of the phase-resolved wave field. A possible consecution of action is as follows: the spectra calculated in the spectral model in the selected areas are used for generation of the initial conditions in the wave number space; then the calculations up to the stationary statistical regime with HWM are performed; finally, the statistical characteristics of wave field including the probability distribution of wave height are calculated. The input and dissipation terms can be included in such simulation. Note that such approach allows calculating nonlinear interactions between all the modes based on full nonlinear equations. Such scheme can significantly extend the output of spectral modeling.
It is comparatively easy to transform the scheme into a finite-difference form and use the model for investigation of a wave regime in a specific geographical environment with variable boundary conditions. Note that a numerical scheme for this model is much closer to the local scheme than FWM, which makes convenient to carry out the parallel computations. An investigation of applicability of a developed approach for the waves above finite depth is underway.