On Quadruplet Interactions for Ocean Surface Waves

Nonlinear wave-wave interactions among ocean surface waves are dominated by quadruplet wave-wave interactions. Computing the nonlinear 4-wave interactions with the Boltzmann integral requires finding the loci of interactions for the quadruplets or solving the delta functions. This is an important part of the computation, but so far it is done by an iteration method that consumes computational time and may not converge after doing several iterations. In this paper, an explicit methodology to find the loci of the quadruplet interactions is presented. This research target is to develop a better method to compute 5 the loci. To illustrate the method, there are 4 cases that will be discussed in this paper. Results show that the new method gives better results than the previous methods that have been applied. Moreover, without iterations the presented method requires less computational loops and some variables, for example the distance between loci, denoted ds, can be determined without any looping. Therefore, the new method leads to better and faster computations than the previous iteration method.


Introduction
In modeling wind generated waves, there are three major factors that must be taken account, which are wind input, dissipation, and non linear wave-wave interactions.Mathematically, these sources are written as, respectively, S IN , S DS , and S N L .See The SWAMP Group (1985), Tolman (2009), Perrie et al. (2013), and Tolman et al. (2014) for details.The first term gives the energy to the waves from the wind, the second term reduces the energy due to wave breaking dissipation, and the last term is conservative, transferring energy among the wavenumbers of the spectrum, Komen et al. (1994), Cavaleri et al. (2007), andHoltuijsen (2007).In general, nonlinear wave-wave interactions (S N L ) among the ocean surface waves are dominated by quadruplet wave-wave interactions, but in shallow water, triad wave-wave interactions may also exist.
Computing the nonlinear 4-wave interactions can be done by so-called exact methods (using the Boltzmann integral) as described by Tracy & Resio (1982), Resio &Perrie (1991), andVan Vledder (2006).Alternately, for operational forecasts we can use approximation methods, such as DIA, the discrete interaction approximation (WAMDIG (1988)).In practical operational marine forecasts, wave models usually use the DIA method, Tolman (2009), but the exact method is also important as a benchmark to solve the problem when the approximation method fails to accurately predict the wave spectra due to unusual wind inputs (for example: rapidly developing storms or hurricanes).
In this paper, we discuss the exact method, partially solving the non-linear interactions with the Boltzmann integral.However, we do not derive all equations of the integral, but focus on a small part of the integration, the important part, the starting point, the domain of the integration.We begin with the main equation, the Boltzmann integral, as shown on Eq. 1.Thereafter, the methodology to solve the equation is explained.Moreover, solving the integral is beyond the scope of this paper.That is done in many other places, for example Resio & Perrie (1991).However, for any definite integral, we have to integrate it over its domain.Regarding the computation, the Boltzmann equation integral is approximated by a collection of nested do-loops; the smaller step the increment in wavenumber space, dk, the better is the result.But on the other hand, when dk is quite small, the computation time becomes excessively large.However, the computation time can be reduced if we do not compute the unnecessary terms (e.g.terms that are zero, or close to zero).The objective of this study is to determine the domain of the integration only, which is the basis to efficiently solve the Boltzmann integral.The content of this paper is to explain the new method mathematically and physically.
As defined by Webb (1978), the Boltzmann integral of the nonlinear action density transfer for quadruplet nonlinear interactions may be written as: where k i are vectors.Thus it is a six-fold integral in terms of k x and k y .It is too costly to use in operational wave forecasting and therefore this term must be approximated if it is to be implemented in practical applications (see Perrie et al. (2013)).
As you see, the integration contains delta functions δ(X).The property of this function says that the integration of δ(X) is zero if X = 0 and one if X = 0, so we only compute Eq. ( 1), on the domain where the following conditions are satisfied Therefore, determining Eq. ( 2) and Eq. ( 3) is the important first part of the process to compute the Boltzmann integral.Points that satisfy those conditions are loci of the quadruplet wave-wave interactions.This is the main topic of this research, finding the loci or domain of the integration.The new method is expected to give better performance than the previous method.
In wave models, usually k 1 and k 3 are given as incident wave numbers, and then we find k 2 and k 4 using the dispersion relation where g is the gravitational acceleration, d is water depth.
In some circumstances, finding k 2 and k 4 is not straightforward.Existing computational methods use iteration schemes, something like where is some error, and the computation is stopped when ≈ 0. This method can take excessive time to solve the dispersion relation Eq. ( 4) and sometimes the iteration does not converge.The computational time can be reduced if we can solve this problem without iterating.Tracy & Resio (1982) shows how to reduce the computational time by using a specific grid geometry.Their work shows that the integrand of Eq. ( 1) is the product of functions: F(action density) and G(geometry term).In polar coordinates, the grid geometry of Tracy & Resio (1982) allows the spectrum to be expressed in terms of wavenumbers k i such that k i = λk i−1 .
Thus, it can be shown that, The computation time is diminished, because the iteration is done for the basic geometry (G 0 ) only.However, the spectral grid is not flexible, because it is a function of the multiplication factor λ. In some applications, it may be desirable to design the grid so that the spectral peak is always exactly at a specific wavenumber, for example the 15th wavenumber (k 15 ), and there are always 15 spectral bins to represent the spectral forward face (k ≤ k 15 ), as the spectrum evolves in time (Resio & Perrie (1991)).However to simulate this case, the basic geometry (G 0 ) needs to be computed at each time step, using the grid geometry of Tracy & Resio (1982) because the spectral peak downshifts to lower frequencies with time.Clearly this methodology is very expensive.

The explicit methods
This is the main section of this paper.Finding loci of the interactions are explained thoroughly.We are to going to solve Eq.
(2) and Eq. ( 3) without iteration to illustrate the method.There are 4 cases that will be discussed.These cases are functions of vector p, angular velocity q, and a dimensionless variable kd.
The defined variables are p and q, where p is a vector, stated as and its components may be written as, and the magnitude is The q variable is an angular velocity difference As you see, the vector p is a function of k 1 and k 3 , so the four cases we are going to consider are determined by vectors k 1 and k 3 : 1. First case, where vector k 1 equals to vector k 3 .
2. Second case, where the magnitude of vector k 1 equals to the magnitude of vector k 3 .
3. Third case, where vector k 1 does not equal to vector k 3 and in deep water.
4. Fourth case, where vector k 1 does not equal to vector k 3 and in shallow water.
Following subsections will show each case thoroughly.Readers may skip the derived equations, they can find the main idea of each case which is summarized at the end of each subsection in italics.Please keep in mind that the derived equations for Case III or IV are built by defining that p is lying on the x-axis as summarized in Section 3.
In this case, the p and q are zero, and k 2 = k 4 .This case will give infinite pairs of k 2 and k 4 , but the term ] will be zero, because n 4 = n 2 and n 3 = n 1 .Therefore, we do not need to find any loci of k 2 and k 4 for this case.

Case II, |k
The variable q will be zero but p = 0.For this case, the loci are not a function of depth, d.We observe that the loci of k 2 and k 4 are normal lines to the slope of p, and some methods bound the normal lines by a circle with radii k max = f (4ω 1 ), to limit the computation.
The slope of p is stated as, then the slope normal to the p will be Therefore, the normal line is defined as, For this case (x 0 , y 0 ) = (k 3x , k 3y ), and so As you see, the loci of this case are lines that are perpendicular to vector p.
In the next section, the loci produce an egg-like shape, as shown on Figure 1.In this paper, the terms 'circle', 'ellipse', and 'egg-like' are interchangeably used to describe the loci of Case III and IV.It is called a circle because it has a center and radii and a starting shape to get the loci.It is called an ellipse because it has a major axis and this shape mimics the real loci.Finally, it is called an egg-like shape because the loci shape really is like an egg.

in deep water
The loci of k 2 and k 4 are on an egg-like shape with its major axis parallel to the vector p.Therefore, we have to find points on the 'end points' of the maximum secant line to get its diameter.These 'end points' are called r min and r max respectively.

Finding r min
In the dispersion relation Eq. ( 4), we have tanh(kd) ≈ 1 when kd ≥ 3, corresponding to the deep water approximation.Therefore, we rewrite Eq. ( 4) in a simpler way From the definition of q in Eq. ( 11) and from Eq. ( 16), we have From the definition of p in Eq. ( 7) and taking its magnitude and the magnitudes of the k i vectors, and remember that these vectors are lying on the x-axis, we can get Therefore, we can write, Let k 2 = r min p, so Eq. ( 19) becomes Equation Eq. ( 20) gives negative values for k 4 when the value for r min is less than 1.However, we need only be concerned about the magnitude of k 4 .Thus we redefine k 4 as We substitute k 2 = r min p and Eq. ( 21) into Eq.( 17), which gives we then take the √ p out of the bracket and do some algebraic re-arrangement, which gives Therefore, we get a quadratic function, r 2 min − r min + C = 0, where We solve for the roots of Eq. ( 23) For the computation, we define r min ≥ 0.5, which implies that we should therefore use the bigger root, which is

Finding r max
For k 2 = r max p, we can write, We substitute k 2 = r max p and Eq. ( 27) into Eq.( 17) which gives and with some algebraic re-arrangement, we find and we simplify Eq. ( 29) as r 2 max − r max = r max − C, where Then squaring both sides of Eq. ( 29), we get

Finding the loci
After r min and r max are found, now we assume that the loci of k 2 lie on a circle-like shape with a center and radius, c c and c r , respectively If we slice the circle-like structure of the loci into n parts, and if the angle of each slice is dφ, then using the triangle identity, we compute k 2 as Therefore, we know the magnitude of k 2 .To find its direction, we must use Eq. ( 17), and rewrite it as, q = √ g( . Therefore, we will get where 2 .We continue with the algebraic derivation by defining the p value which is lying on our x-axis, (p, 0).Therefore, we obtain The above equation specifies another circle with the center at p and radius ζ.We continue with the solution of Eq. ( 36) and substitute k 2 in polar coordinate form; (k 2 cos θ, k 2 sin θ).Therefore, we get Therefore, the angle between k 2 and our axis is defined as Finally, using Eq. ( 34) and Eq. ( 38), we can determine the magnitude and direction of k 2 , and once k 2 is determined, the k 4 is found.Loci of the tips of these vectors depict an egg-like shape.
As one can see, in this case we assume tanh(kd) ≈ 1, to simplify the problem.We compute the minimum point (tip of the smallest vector) and maximum point (tip of the largest vector).Next we find the center of this interval, then other points (tip of vectors in between) surrounding the center can be determined.

Case IV, k 1 = k 3 , in shallow water
As in Case III, the loci of k 2 and k 4 are also an egg-like shape with its center parallel to the vector p.However, the points on the 'ends' of the maximum secant line are found by a linear approximation method, in this case.

Finding r min
Because the kd is not large compared to 1, therefore we cannot neglect the term, tanh(kd).Thus, from Eq. ( 11) and Eq. ( 4) and substituting values, k 2 = r min p and k 4 = (1 − r min )p respectively, we get The square root is a multiple products function.To find r min is not a straightforward problem.Previous efforts tend to try to solve the problem without simplifying it first.Most of these approaches use an iteration method to find the unknown roots.An example is the computation done by Resio & Perrie (2008).
In this study, we use another approach, which is straightforward and without iteration.In this approach, we propose an explicit analytic solution, which can find the solution easily and much faster than the previous iterative approaches.
Let's go back to Eq. ( 39).Before we solve this equation, we focus on the dependent variable, r min .
Let's call this variable x.Therefore, we get We rewrite this equation in a simpler way where H = f g.Now, we can see clearly what we have.The q variable is just a difference between 2 functions.In fact they are the same function, but the second term on the right hand side is a mirror, f (−x), of the first one, f (x), and also shifted +1, as f (1 − x).Let's examine the behaviors of H(x) and H(1 − x) and we will plot them in the domain 0 ≤ x ≤ 1.We find that both equations are linear, and therefore we can conclude that the difference is also linear.For example, if we have linear functions, where m c = m 1 − (−m 1 ) = 2m 1 .The difference is zero, when r min = x = 0.5, which means that the function is shifted to the right.Therefore the difference function, q, may be stated as Replacing x with r min , we get This method is depicted on Figure 3.We find that the difference of two functions, q, is linear, and so we have to find its slope and shifting factors.Once we have the equation q = f (x), then the variable x (or r min ) will be known, because q is provided.
Figure 3. Finding the rmin.The q = f (x) − f (1 − x) behaves linearly on domain x < 1.Therefore, if q is given, then x or rmin can be found.
-Based on these inputs, the program will select the case, from the four described above.
-For Cases III and IV, the computation is done on the x-axis first (assume p(p, 0)), and then the solved loci are rotated by θ p degrees anticlockwise, where θ p is in the direction of p. Figure 5 depicts this mechanism.
-Finally the loci of k 2 , k 4 , and ds are determined, where ds is the distance between two consecutive points.The derived equations on the previous sections are based on the assumption that k 3 > k 1 .If k 3 < k 1 , the integral is the 10 negative of the integration that is computed by the k 3 > k 1 condition.
The previous method by Resio & Perrie (2008) uses 50 calculations to iterate the r min and an error flag if the computation does not converge.To find the r max , it requires 200 calculations for a coarse approximation and 10 additional calculations for refining it, and it is also equipped with an error flag.However, the new method requires no iteration and converges.Thus the computational time is faster than that of the old one.
The results shown here are for the loci of k 2 and the resonance condition, another disadvantage from the iteration method, besides being slow.Figure 7 shows results of both methods on the same axis for comparison.In the Case II, where p = 0 and q = 0, we should get a normal line to the vector p.Although the difference is small, the old iterative method cannot produce a straight line as expected, especially when p is horizontal or vertical.See Figures 7a and 7b for details.
For other cases, the iteration method gives a small residue (about O(10 −8 )) for the sum of k 1 + k 2 − k 3 − k 4 in Eq. ( 2), but the new method can achieve this condition more accurately, as shown in Figures 8a and 8b, essentially zero.Now, we discuss the result of the requirement that the frequencies sum up to a total of zero (Eq.( 3)).As explained in Section 2, finding that the summation of the wave numbers should equal zero involves q (Eq.( 11)), the angular velocity difference.
Therefore, once the loci are found, we can conclude that the total sum of angular velocities is also zero.and ω 2 = ω 4 regardless the water is finite depth with tanh(kd) or not.A sample of resonance conditions for Cases III and IV are tabulated in Table 1 and Table 2 respectively.
Because the solution is straightforward and because no iteration is needed, setting up the grid is not necessary, because the new method can work for any grid.This is clearly an improvement over the iterative method.The iteration method (labeled 'old') produces small a residue for the sum of kn, but the explicit method (labeled 'new') gives no residue.Only sum of vector components is presented here.

Conclusions
It is shown that both methods can determine the resonance conditions, k However, the new method gives less residue.The residue is a number (its magnitude) of adding all solved wave numbers, this number should be zero.Therefore, the smaller the residue the better the method.The new method also omits iteration so it   On the other hand, the explicit method can simplify the computation.For example in the Case I, to determine the distance between locus, ds, the old method will wait for all loci, locus 1 , locus 2 , . .., to be first found.Then it will do subtractions, ds 1 = locus 2 − locus 1 , ds 2 = locus 3 − locus 2 , etc.However, the new method can determine the ds at the beginning without waiting until all loci computed.For example, because we know that in Case I, the loci is a straight line, then the ds i = length of the line / (number of locus -1), and ds is uniform.
The existing wave model, for example SWAN or WAVEWATCH III T M , computes the S N L based on the deep water (Case III) formulation.Applications where water is shallow of finite depth need to check again how the computation of S N L is derived.In this paper, we show 4 cases to determine the loci or domain of integration.Four subroutines to compute the loci can be easily built to handle any condition (deep or shallow water).Next, after we get the loci correctly, we have to examine other variables in Eq. 1. Are the subprograms derived from formulations that suit the condition?However, actually finding the other variables is beyond this paper.This paper only shows how we start to do the computation for S N L using full Boltzmann integral correctly, by finding the loci based on the water depth.A correct procedure that is followed at the beginning of simulation will produce good results.
If we have proper subroutines that suit any condition of the real sea, then we can do simulation for hindcasts or forecasts with better results.The applications to societal concerns and problems are numerous.Good simulations and forecasts of ocean wave data can make good contributions to studies related to ocean wave energy resources utilization, disaster prevention and reduction, and search and rescue, etc (Kamranzad et al. (2013); Mirzaei et al. (2013); Zheng et al. (2014); and Zheng & Li (2015)).

Code availability
The code written in FORTRAN is available on the following link https://doi.org/10.5281/zenodo.1040466

Figure 1 .
Figure 1.An example of loci from Case III or IV.
Figure 2.After the minimum and maximum vector are found, other vectors in between are determined.

Figure 5 .
Figure 5.The sketch of the mechanism to get the loci.

Figure 7 .
Figure 6.Results from different methods for Case II with p have the some slope.The iteration method, graph (a), gives non uniform ds and exceeds the boundary 4f .On the other side, graph (b), the explicit method gives uniform ds and does not exceed the boundary.

Table 1 .
Example of the resonance of ω in the deep water case.