Clarification of the 2-D Formulation for 3-D Potential Periodic Waves

The advantages and incompleteness of 3-D direct (phase-resolving) wind wave model (named as TriDWave model) are discussed. The refinement of a two-dimensional formulation for the three-dimensional phase-resolving model of potential waves is proposed. Both models are used to simulate development of the wave field with the same initial conditions over one hundred thousand time steps. The applicability of a 2-D approach is confirmed by comparing the evolution of some integral characteristics and the probability of the elevation. An accelerated model runs approximately 12-15 times faster than the full one. Both models are intended for studying statistical quasi-stationary or non-stationary regime of multimode wave field. The full model can be used for reproduction of periodic wave processes of any complexity, when a single-valued surface is preserved.

Keywords:Phase-resolving modeling; 2-D formulation for a 3-D problem; Evolution of integral characteristics; Wave height probability

The numerical modeling of periodic adiabatic waves at constant depth is not a difficult task. In practice, such a regime can be reproduced for a relatively short time or artificially maintained by including small energy input to compensate the high-frequency damping introduced to provide stability. This version of the model can be used to study the statistical characteristics of a wave field for any initial conditions that ensure stable calculations over a sufficiently long period. Naturally, the statistics should be obtained after the model has passed relaxation and entered an equilibrium regime. A much more difficult task is to reproduce the evolution of a wave field in the presence of energy input and dissipation [1]. Unfortunately, this part of wave theory, despite the abundance of scattered data and a large number of publications, still is full of contradictions and does not provide clear recommendations. The most obvious type of dissipation is filtering of wave disturbances in the spectrum tail [2]. Several options can be proposed here, but the goal remains the same: reducing the accumulation of energy to prevent computational instability which develops along the physical instability. The rate of such dissipation should be minimal, so that the tail does not absorb energy from the main part of the spectrum.

It follows that development of methods for parameterization of physical processes in waves is far from being complete. Following the overturning of waves in the sea, the waves survive, although in a modified form. A numerical model always stops without the possibility of resuming the calculation after the overturning. An algorithm describing wave breaking should start working shortly before the overturning and ensure further calculation at a reasonable level of energy loss. A long-term evolution of wave spectra is completely determined by the spectral distribution of the combine effect of input and dissipation of energy (1). The experimental data does not allow us to determine the role of those processes separately. For example, we can consider the phenomenon of wave breaking which normally occurs at wave peaks. Since it is a discrete process, its spectrum will certainly be shifted to higher wave numbers, probably, close to the region of high-frequency damping. In fact, it is the height of the nonlinear mode that decreases. It follows that the parameterization of breaking cannot be formulated in Fourier space. The process develops in a physical space where the object is the surface itself, not its distorted image in Fourier space. In principle, a detailed study of the overturning process can be formulated and investigated in terms of a two-phase fluid, but such approach is certainly not applicable to long-range modeling of an evolution of a multimode wave field.

The most interesting task which can be solved with a 3D phase resolving model is the development of a wave field under the action of wind [1]. This process is accompanied by the input of energy and dissipation distributed unevenly and overlapping with each other in the spectrum. The growth of the integral wave energy reflects a combined effect of those processes. The primary task of modeling is to accurately reproduce the change of the integral wave energy as the wave field develops. This task, as applied to the problem of wave forecasting based on the spectral models, is solved by choosing suitable parameters based on broad comparison with the observational data. Since the spectral models consider just the density of wave energy as a function of direction and frequency, the problem of separate determination of the role of energy influx and dissipation is usually not posed. The total effect of those processes is actually considered. With a fixed scheme used to calculate energy input, the parameters determining the spectral distribution of a dissipation rate are to be tuned. The phase-resolving models opposite to the spectral models reproduce the evolution of primary variables, i.e. the velocity field (velocity potential) and the shape of the surface. In fact, in such models there is no direct need to use Fourier representation of variables, since they are treated as the continuous functions specified in a grid form. In the periodic domain, Fourier representation is used to improve the accuracy of differentiation and to save the number of operations, since the dimensions of an array of Fourier coefficients are usually several times smaller than those for the array of the grid values. The advantages of linear operations in Fourier space over those in a grid space are obvious. For example, the nonlinear interactions between Fourier modes are usually calculated on a limited set of modes, while in the grid space all the interactions allowed by the resolution are taken into account. The Fourier representation is, of course, convenient when processing the results.

Developing the potential and adiabatic versions of the 3-D direct model is a relatively simple problem of computational mathematics. Such a model can be used for the short-period calculations to calculate the statistical characteristics of waves with a fixed wave spectrum. However, such a model cannot be considered as a model that reflects real properties of wind waves which are formed and evolve solely under the influence of energy input and dissipation. It is known that developing algorithms for calculation of those processes requires multiple repetitions of runs in order to select a specific scheme and parameters. A full model includes a 3-D equation for the velocity potential, the solution of which takes 95% of the time. Two years ago we encountered the practical impossibility of working on the physical blocks of the model due to excessive expenditure of human time. The peculiarity of a wave model is that it does not work properly when a low spectral and vertical resolution is used. Our experience shows that the number of degrees of freedom for surface variables (the number of Fourier coefficients for the level and potential velocity) cannot be significantly less than 106 . The decrease in vertical resolution leads to rapid increase in the error when calculating the vertical velocity on the surface. Such difficulties helped us to realize the necessity of developing an accelerated version of a model capable of generating the results close enough to those obtained with a full model. The idea is based on the use of a 2-D equation for the vertical velocity instead of a 3-D equation for the velocity potential. Several versions of an accelerated model together with demonstration of the verification results are presented in the works [3-6]. The most successful scheme is described in this article. The scheme seems suitable for developing some methods for parameterization of the physical processes in waves; therefore, its principal modification is not planned in the nearest future.

A new algorithm is based on projection of an equation for the 3-D velocity potential onto the surface ζ = 0 (see (4)):

All variables here are nondimensional on the basis of the length scale L and the gravity acceleration g. The representation used here is w = w+ w where w is a full vertical velocity; w is its linear component calculated by a formula

where τ = t and ξ ,ϑ,ζ are the new coordinates connected with the Cartesian coordinates by relations

η η are the derivatives of η , w is a linear component of vertical velocity, which Fourier amplitudes are calculated with a formula

are the Fourier transform basic functions (see ( )),

In the articles (3-5 ) it is shown that the vertical derivative of a nonlinear component w_ζ  in (1) is related by a linear relation with the nonlinear component w itself:

where σ_η is dispersion of the elevation η while the coefficient A is a function of an additional, truly nondimensional (i.e. not depending on external length scale L) parameter μ

The dependence A(μ ) was investigated on the basis of a full 3-D model TriDWave. Firstly, the evaluation of this function was attempted in a Fourier space [3]. Such a method turned out to be unsuccessful because of large scatter of values. That result was expected since we have seen many times that the representation of a field as a collection of a large number of linear modes with random phases turns out to be inconsistent with reality. It works only for a traditional assumption of small nonlinearity of wave field. In reality the physical modes are far from being linear, and they change their properties quite quickly. The primary continuous fields approximated by the grid variables are formed by a large number of modes, their statistical properties being much more stable than those for Fourier amplitudes. In general, direct models do not consider the amplitudes of modes as the physical variables at all. The Fourier representation is only an element of a Fourier transform method and a traditional way of processing the results. For wave modeling in an arbitrary domain it is more convenient to consider only the grid variables.

For evaluation of function the results of modeling of development with a TriDWave, having a grid resolution 1024×512 and a number of levels equal to 20 were used. A description of a similar formulation of the problem is given in (4). During the experiment the pairs (w,w_ζ) of values were recorded. The instantaneous fields containing 524,288 of such pairs were used to determine the coefficient A in the linear dependence [5]. The number of such records was equal to 99, hence the total number of pairs was 51,904,512. The dependence of the obtained values of the coefficient on the parameter A is shown in Figure 1 by a clearly defined grey area.

Figure 1:Dependence of A on μ =σ_η σ_L . The data calculated with a full model are indicated by grey dots. A solid curve corresponds to an approximation (6) extended above and below the area of the calculated values c A . The calculated data closely coincides with the approximation; the mean difference between them is about 0.03.

As seen, the data A lies on a smooth curve with no deviations. Increasing the parameter μ means increasing the nonlinearity. The approximation of this dependence is

An equation for calculation of a non-linear component of the vertical velocity takes the following form

Note that the right-hand side of the equation contains w since w = w+ w . The equation is solved by plain iterations until the increment ì àõ δ w exceeds ε . In the present calculations he value ε =10−7 was used. The number of iterations normally does not exceed 2; at very large steepness it reaches 4. А non-convergence of iterations as a rule always indicates a non-single valued surface occurrence.

The reliability of a 2-D model can be most convincingly demonstrated by the point-to-point comparison of the elevation field. A correlation coefficient C and the normalized root-meansquared difference between the fields

The coincidence of the elevation fields decreases with time, however, over a time of about 200 periods of the peak wave (what corresponds to 13,000 time steps), the correlation between them dropped to just 0.95. Considering that the surfaces with the number of degrees of freedom of the order of 10⁶ have a complicated multimode quasi-random nature, it was expected that violation of the similarity would occur a lot faster. We consider the result shown in Figure 2 as quite satisfactory.

Figure 2:Comparison of the evolution of the elevation fields simulated by the full 3-D and accelerated 2-D models: T_p is the time expressed in periods of an initial peak wave; C is a correlation coefficient; D is a normalized root-mean squared difference (Eq. 8).

The formulation of the 3-D phase-resolving models based on full equations requires a detailed description of energy transformation. As mentioned above, these problems are not sufficiently investigated. The author believes that some progress has been made in the physics of the wind wave dynamic interaction. The description of energy influx to waves is based on Miles (1957) theory [7] which states that the complex amplitudes of each mode of the surface are connected to those of the surface pressure by a linear relation with a complex β -coefficient approximated as a function of some nondimensional frequency. Miles theory initially formulated for a wind field consisting of the small amplitude modes turns out to be surprisingly flexible. The real field consists of superposition of Stokes-like waves characterized by numerous nonlinear features such as sharp crests, gentle troughs, group structure and others. The air flow above them contains multiple flow separations, shadow zones and other transient patterns. The function β underlying the parameterization of the energy input to waves was found as a result of massive calculations with a model describing the dynamic interaction of waves and wind [8]. In a cited work a reliable scheme was formulated. Qualitatively, that scheme agrees with early version used in WAVEWATCH model [9], but it predicts slower increase β -function with a frequency. It is suitable for including in spectral prognostic and phase-resolving models. The scheme for calculation of input energy to waves is to be further developed. However, the situation with the wave breaking problem remains unresolved. A developed 3D model based on full equations in principle is able to simulate development of the wave field characterized by growth of energy and downshifting. Such numerical experiments were carried out and the results seemed quite plausible. However, in each case it was unclear how realistically the energy growth and downshifting rate were reproduced. For example the JONSWAP approximation [10] predicts that growth of energy is proportional to τ^3/2 . This rule is valid probably to very short time τ . As estimates and calculations give, the initial rapid growth of energy under the influence of increasing dissipation is replaced by a much slower rate and eventually settles at the level estimated by the Pierson/Moscowitz spectrum [11].

The modern version of the wave breaking energy dissipation algorithm is based on highly selective smoothing of the elevation field and surface potential. Under the constant wind action, the wave field structure becomes more complicated: The energy, steepness and curvature of the surface increase. Eventually, instability occurs. The survival period can be extended by decreasing the time step, but only for a short period. If a non-single-valued surface occurs somewhere in the region, the model terminates immediately. If the task is to simulate a further wave development, the instability must be prevented by smoothing. This algorithm has been repeatedly tested and described (see [2]). The evolution of integral characteristics of the solution during the first 100,000 steps is shown in Figure 3.

Figure 3:The evolution of integral characteristics calculated over 100,000 time steps. A straight line 1 indicates the zero value of integral transformation of energy due to nonlinear interactions (it is smaller by 5 decimal orders than the rates of energy transformations). 2 – the rate of energy input from wind; 3 – the rate of dissipation; 4 – energy balance (energy input minus dissipation); 5 – total energy. Values 2, 3 and 4 are multiplied by 105, the energy 5 is multiplied by 10⁵. Black curves refer to a full model while red curves – to the accelerated model.

The comparison of integral characteristics proved that the results obtained with different models are close to each other. Thus, the task that was set for construction of a simplified scheme is achieved: an accelerated model runs much faster than the full one and it gives quite acceptable results. In this case, the time spent on the calculations is 13 smaller than that for full model. Figure 3 shows, however, that there are small discrepancies between the integral characteristics. Interestingly, the overestimated energy influx (curve 2) is partially compensated by increase in dissipation (curve 3), so the energy balance (curve 4) is reproduced with a smaller error, which has little effect on the change in energy (curve 5). In any case, this divergence deserves to be studied and corrected. When discussing the results, in the first place it is important to remember the purpose for which direct models are developed. We believe that the main task is to study the statistical properties of the real wave field: The probability distribution, the spectra of various variables and their dependence on the degree of wave development. In order to successfully solve this problem, it is necessary to refine the parameterization of physical blocks of the model. Figure 4 shows a comparison of the most frequently considered characteristics, i.e. the distribution of elevation probability. The figure demonstrates very close agreement between the results obtained with totally different approaches. Since the curves coincide, the absolute difference between them is shown by dots. Note that the calculations with a full 3-D model took about 12 days, while the same with a 2-D model took less than 1 day.

Figure 4:The probability distribution for the elevation fields calculated over the entire period of integration. A thin solid line corresponds to a full model, while a dashed line – to the accelerated model. These curves nearly coincide. The absolute difference between the data is shown by dots. As seen, the difference is about 3 decimal orders smaller than the probability itself.

The article describes the latest version of reduction of the 3-D model to the 2-D model for 3-D waves. The main advantage of the simplified model is a quite high calculation speed while maintaining the accuracy of reproduction of the statistical characteristics of wave field. In fact, direct wave models are created to solve this problem. We do not propose to abandon direct three-dimensional modeling; on the contrary, the three-dimensional model [2] was created specifically for this work. The problem is that all existing models are not actually the models of the real process. Creating an accurate model for adiabatic waves is a simple task. Its only drawback is that the adiabatic motion does not exist in nature. The structure of the wave field and the nature of its evolution are completely determined by energy input and dissipation. Those processes are still poorly understood. In fact, they are much more complicated than the entire theory of adiabatic waves. The work on development of parameterization of the physical processes cannot be done without multiple long-range calculations with a direct model. Currently, these calculations take too much of computer time for such activity to be successful. A simplified model solves this problem. Meanwhile, due to the fact that the simplified model allows us to reproduce statistical properties of waves with satisfactory accuracy, it can be used for these purposes quite successfully too. The full and accelerated components are combined in one model called TriDWave model. The choice of the scheme can be made by changing one parameter. The paper with a detailed description of TriDWave is in preparation.

Author would like to thank Mrs. O. Chalikova for her assistance in preparation of the manuscript.

1. Chalikov D (2018) Numerical modeling of surface wave development under the action of wind. Ocean Sci 14: 453-470.
2. Chalikov D (2016) Numerical modeling of sea waves. Earth and Environmental Science, Earth and Environmental Science (R0). 1^st (Edn.), Springer Publishers, USA, pp. 1-330.
3. Chalikov D (2023) The 2-D heuristic wave model for 3-D waves with a local closure scheme. Examines Mar Biol Oceanogr 5(2): 1-5.
4. Chalikov D (2022b) A 2D for 3D Periodic Deep-Water Waves. J Mar Sci Eng 10: 410.
5. Chalikov D (2021) Accelerated reproduction of 2-d periodic waves. Ocean dynamics 71: 309-322.
6. Chalikov D (2021) A two-dimensional approach to the three-dimensional phase resolving wave modeling. Examines Mar Biol Oceanogr 4(1): 1-4.
7. Miles JW (1957) On the generation of surface waves by shear flows. Journal of Fluid Mechanics 3: 185-204.
8. Chalikov D (2025) Numerical investigation of energy exchange between wind and waves. Water Waves.
9. Tolman H, Chalikov D (1996) On the source terms in a third generation wind wave model. Journ Phys Oceanogr 11.
10. Hasselmann K, Dirk O (1973) Measurements of wind-wave growth and swell decay during the Joint Sea Wave Project (JONSWAP). Tsch Hydrogh Z Suppl A8(12): 1-95
11. Pierson WJ, Moscowitz L (1964) A proposed spectral form for fully developed wind seas based on the similarity theory of S. A. Kitaigorodskii. Journal of Geophysical Research 69(24): 5181–5190.