Prediction of height and period joint distributions for stochastic ocean waves

This article proposes a new methodology to predict the wave height and period joint distributions by utilizing a transformed linear simulation method. The proposed transformed linear simulation method is based on a Hermite transformation model where the transformation is chosen to be a monotonic cubic polynomial, calibrated such that the first four moments of the transformed model match the moments of the true process. The proposed new approach is applied for calculating the wave height and period joint distributions of a sea state with the surface elevation data measured at an offshore site, and its accuracy and efficiency are favorably validated by using comparisons with the results from an empirical joint distribution model, from a linear simulation model and from a second-order nonlinear simulation model.


Introduction
In the practice of ocean engineering, it is very important to study the statistics of the wave heights and the associated wave periods. Firstly, wave breaking takes place when the wave height and period cannot maintain the equilibrium condition needed for stability. Therefore, it is necessary to have the knowledge of the joint distribution of the wave heights and the associated wave periods in order to predict the possibility of the occurrence of wave breaking in a given sea state. Further, in designing a floating marine system, it is extremely important to know the statistical information concerning the wave height with periods that are close to the system's natural period in areas where the structure will be operated (Ochi, 2005). Furthermore, studying the statistics of the wave heights and the associated wave periods is also very important in wave power studies (Pontes, 1998;Gonçalves et al., 2014;Zodiatis et al., 2014).
In a real world ocean engineering project, the known information regarding the environmental conditions is typically a wave spectrum corresponding to a short-term, stationary sea state (Nazir et al., 2008;Vanem, 2011). Based on this specific wave spectrum, there are several approaches to predict the wave height and period joint distributions. The first and most straightforward way is to resort to some empirical (or theoretical) models, such as those proposed by Longuet-Higgins(1975, 1983, Cavanié et al. (1976), Lindgren and Rychlik (1982), Stansell et al. (2004) and Zhang and Guedes Soares (2016). The joint distribution model proposed by Longuet-Higgins (1975) has been developed based on the narrowband linear Gaussian wave process theory and suffers the drawback that it is symmetrical with respect to the normalized wave period, while the asymmetric property is detected after a detailed analysis of the experimental data accomplished by Chakrabarti and Cooley (1977) and also by Goda (1978). The asymmetric model for the joint distribution developed by Cavanié et al. (1976) is based on the similar assumption of a narrowband linear Gaussian wave process. The analytical joint distribution model for waves with arbitrary spectral bandwidth proposed by Lindgren and Rychlik (1982) is also based on a linear Gaussian wave process theory and the distribution is not given in a closed form, and the computation is too complex to be used in practice. Meanwhile, a common problem of the latter two theoretical models is their critical dependence on the forth-order moment which is strongly affected by the tail of the wave spectral function. As a consequence, computational factors such as the sampling rate and the high-frequency cut-off used to estimate the spectral moments can affect them. Later, Longuet-Higgins (1983) revised his theoretical model from Longuet-Higgins (1975) and presented a new distribution that is asymmetric to normalized period and depends only on the lower-order spectral moments. Stansell et al. (2004) and Zhang and Guedes Soares (2016) made further improvement of the Longuet-Higgins (1983) theoretical model for Gaussian process. However, it can be obviously noticed that all the aforementioned six theoretical joint distribution models have been derived based on the linear Gaussian wave process theory (e.g. Rydén, 2006), which contradicts the fact that ocean waves in the real world are nonlinear and non-Gaussian. Moreover, the most recent studies have also shown that the results from the joint distribution model of Longuet-Higgins (1983) deviate a lot from those obtained from the measured sea data (Shu et al., 2012;Guo et al., 2014Guo et al., , 2016 or from the generated laboratory wave data (Zhang et al., 2013).
From the specific wave spectrum, there is another way to obtain the wave height and period joint distributions by using the numerical simulation method. Based on some idealized target spectra, Rodríguez and Soares (1999) and  performed linear simulations to obtain long wave time series which satisfied the Gaussianity and stationarity hypothesis. Their simulations involved the superposition of a finite number of harmonics with different amplitudes, frequencies and a random phase (Shinozuka and Wai, 1979;Hudspeth and Borgman, 1979). They next aggregated all the simulated waves into a joint histogram to obtain the wave height and period joint distributions. However, their simulated waves (and the subsequently obtained wave height and period joint distributions) are good approximations to the real ones only under the condition of small oscillation amplitudes. Meanwhile, Rodríguez and Soares (1999) and  did not verify their research results with field data. Wang and Xia (2012) used a nonlinear simulation method to compute the wave surface elevations for the second-order irregular Stokes waves. In their numerical implementation, they divided the problem into two processes, i.e. first generating a first order linear (Gaussian) wave process with a specific wave spectrum, then adding the nonlinear corrections to it based on the second-order hydrodynamics. The obtained time histories of waves were then statistically processed to acquire the wave height and period joint distribution. They had utilized the wave data measured at the Gulfaks C platform in the North Sea to verify that their nonlinearly simulated joint distribution is better than that calculated by using the Cavanie et al.'s joint distribution model. However, the nonlinear simulations performed by Wang and Xia (2012) are too time consuming, and this drawback will certainly affect the applicability of the nonlinear simulation method in some time-constrained offshore engineering projects.
In this paper, a new methodology to predict the wave height and period joint distributions by utilizing a transformed linear simulation method will be proposed. The proposed transformed linear simulation method will be based on a Hermite transformation model where the transformation is chosen to be a monotonic cubic polynomial, calibrated such that the first four moments of the transformed model match the moments of the true process. The proposed new approach will be applied for calculating the wave height and period joint distributions of a sea state with the surface elevation data measured at an offshore site, and its accuracy and efficiency will be verified by using comparisons with the results from an existing joint distribution model, from a linear simulation model and from a second-order nonlinear simulation model.

Nonlinear Stokes waves
In an idealized Gaussian random sea model the individual cosine wave trains superimpose linearly (add) without interaction. Therefore, the model is also called the linear sea model. However, it is known that waves in the real sea are nonlinear, and the ocean surface elevation process deviates from the Gaussian assumption, i.e. the wave crests will become steeper and higher and the wave troughs will become flatter and shallower than expected under the Gaussian assumption (Ochi, 2005). In the following we briefly introduce the theories of the nonlinear Stokes waves.
Φ (x, y, z, t) Φ (x, y, z, t) For an idealized flow, the flow velocity can be described as the gradient of a velocity potential in which x and y are the horizontal Cartesian coordinates, z is the vertical Cartesian coordinate and t is time. The positive z-direction points upward. The free surface is located at z=η (x, y, t), and the bottom of the fluid region is at z=-d(x, y). Assuming incompressible flow, for constant water depth d, the potential and the surface elevation η (x, y, t) are determined by the following boundary value problem (e.g., Toffoli et al., 2006): (1) With a perturbation-series approach (Stokes expansion), the solutions of the system Eqs. (1)-(4) can be sought using the following equations (Toffoli et al., 2006): In Eq. (5) is a small parameter in the expansion and it is typically proportional to the wave steepness. For an irregular sea state characterized by a specific wave spectrum where is the angular frequency, it is straightforward to show that a first order linear Gaussian solution of the system Eqs.
(1)-(5) takes the following form (Toffoli et al., 2006): as N tends to infinity. In Eqs. (6) and (7), Re denotes the real part of the complex number, and for each elementary cosine wave c n denotes its complex valued amplitude which is directly related to the specific wave spectrum by the relation . In Eqs. (6) and (7), is the angular frequency and , and is the upper cut off frequency beyond which the power spectral density function.
may be assumed to be zero for either mathematical or physical reasons. is the wave number, and is the phase angle that is uniformly distributed in π ω n k n the range of [0, 2 ]. Meanwhile, in Eqs. (6) and (7) the individual frequencies, and wave numbers, are functionally related through the so-called linear dispersion relation: where g and d are the gravitational acceleration and water depth, respectively. By numerically implementing Eqs. (7) and (8), we can generate the first order, linear (Gaussian) time histories of the sea waves at a specific location, and this technique is usually called a linear simulation method in the existing ocean engineering literature.
Waves in the real ocean usually do not follow the linear Gaussian model. The linear Gaussian sea model can be corrected by including the quadratic (second-order) terms. Following Langley (1987) and Hasselmann (1962), the quadratic (second-order) corrections are given by Eqs. (9) and (10), where , r mn and q mn are called the quadratic (second-order) transfer functions. The quadratic transfer function is given by Eq. (11), in which the Kroenecker delta ( =1 if n+m=0, zero otherwise) is introduced to avoid a singular . The quadratic (secondorder) transfer functions r mn and q mn are given by Eqs. (12) and (13). ; (11) ω n In Eqs. (12) and (13), the wave numbers k n and frequencies satisfy the same relation as in the linear case. Fi-nally, by combining Eqs. (7) and (10), the wave surface elevations for the nonlinear Stokes waves can be written as: By numerically implementing Eq. (14), we can obtain the wave surface elevations for the second order irregular Stokes waves, and this technique is usually called a nonlinear simulation method in the existing ocean engineering literature. For the numerical implementation of the nonlinear simulation, the problem is usually divided into two processes, i.e. first generating the first order, linear (Gaussian) time histories with a specific wave spectrum at a specific location, for each of these numerical algorithms then evaluates the full set of the second-order corrections according to the hydrodynamic theory as specified in this section. Therefore, the first order wave process, with N components at frequencies , gives a rise to a total of N 2 corrections, spreading over all sum frequencies at frequencies , and to another N 2 corrections over all different frequencies . However, the nonlinear simulation method as explained above will become very time consuming when N becomes large. In order to overcome this drawback, we propose to utilize a transformed linear simulation method in this article to predict the wave height and period joint distributions, and in the next section the theoretical background of this method will be elucidated.

Introduction on the transformed linear simulation method
An alternative and faster method rather than including the second order correction terms to the linear model, will be implemented by using a deterministic and memoryless functional transformation. The non-Gaussian process, is then a function of a single Gaussian process, : where G(·) is a continuously differentiable function with a positive derivative. In this paper, the function G(·) is proposed to be based on a Hermite transformation model. This transformation is chosen to be a monotonic cubic polynomial, which is calibrated such that the first four moments of the transformed model match the moments of the true process. In the following, we will write as . Let us first give the definition of the Hermite polynomials. The Hermite polynomials of degree n, denoted by , is defined as a function to satisfy the relationship given by We next define the notion of a "Hermite moment" h n as (Wang, 2014): (17) We can then use the following formula to calculate the function G(·) (Wang, 2014): In the above equation, and are respectively the mean and standard deviations of the surface elevation process . The coefficient is a scaling factor, while the coefficients may be related to the Hermite moments in Eq. (17) by applying a Hermite polynomial to Eq. (18) and take expectations. For N=4 moments, the solution is (Wang, 2014): In the above two equations, and are respectively the skewness and kurtosis of the surface elevation process . Eqs. (18)-(20) are the core procedures within the "transformed linear simulation method" in this paper. Meanwhile, we can obviously see that the accuracy of the functional transformation calculated according to Eq. (18) depends on how many moments are kept during the calculation. The more moments are kept, the more accurate the functional transformation will be.
From the above mentioned derivation process it can be noted that in order to calculate the function G(·), we need to know the values of the skewness, kurtosis, mean and standard deviation of the second order nonlinear random waves. In the following, the detailed mathematical procedures to obtain these statistical parameters will be illustrated.
Without loss of generality we set the position coordinate x in to be x=0. Meanwhile, we also write in a simplified way as to be . Then, after a lengthy derivation process as shown in Langley (1987), Eq. (14) in Section 2 can be re-written in the following matrix notation: where Q and R are the real symmetric matrices whose nm-th components are and , respectively; s is a vector whose n-th components are s n ( ). X and Y are vectors whose n-th components are x n and y n , respectively; and x n and y n are the Gaussian random variables whose expressions are shown in Langley (1987). By performing an eigenvalue decomposition, the above equation becomes: where is a diagonal matrix with the eigenvalues of the matrix in the respective diagonal and P 1 contains the corresponding eigenvectors per row. Meanwhile, is a diagonal matrix with the eigenvalues of the matrix in the respective diagonal, and P 2 contains the corresponding eigenvectors per row. By introducing a new set of Gaussian random variables Z j , such that: η(t) we can then write the stochastic process as (this is called the Kac-Siegert solution): where Z j are the independent Gaussian processes with unit variance, and and are the coefficients computed based on the information provided by the sea spectrum , which is chosen for a given sea state. Specifically, the coefficients in the above equation are computed as follows: where denotes a column vector formed by the diagonal elements of . Having obtained the values of and , the values of the mean, standard deviation, skewness and kurtosis of the second order nonlinear random waves can then be calculated as (Langley, 1987):

Calculation examples and discussions
In this Section we will utilize our proposed transformed linear simulation method to calculate the wave height and period joint distributions of a sea state with the surface elevation data measured at the Poseidon platform in the Japan Sea. Poseidon platform was located 3 km off the coast of Yura in the Yamagata prefecture in the Japan Sea during the measurements (Yago et al., 1991). Fig. 1 shows the geographical position of the coast of Yura. The surface elevation data at the Poseidon platform were measured with three ultrasonic wave gauges located at the sea floor from 8:12 a.m. on Nov. 24, 1987 to 7:57 a. m. on Nov. 25, 1987 and three sets of data each consisting of 85547 wave elevation points had been obtained during the measurement. The sampling rate is 1 Hz. The water depth at the Poseidon platform is 41 m. The area was chosen due to a severe seasonal wind from West-North-West during the winter period. In our research, we have divided each of the above three sets of data into eight short records respectively when we calculate the wave height and period joint distributions. In this study we first choose the first three-hour record of the wave elevation points time series measured with #2 ultrasonic wave gauge. Please note that these data are used for academic research only but not for commercial purposes. In Fig. 2, the power spectrum corresponding to this specific three-hour wave record is shown. Fig. 3 shows the calculated wave crest height and period joint distributions for this measured specific three-hour   sea state. The solid black lines in Fig. 3 represent the results of the wave crest height and period joint distributions directly obtained from the measured three-hour wave elevation points at the Poseidon platform. To obtain these results, the time series of the wave crest height are first extracted from these three-hour wave elevation points. The time series of the wave crest front period and crest back period are also extracted from these three-hour wave elevation points, and adding these two together we can obtain the wave period time series. Then exact Epanechnikov kernel density estimates are carried out to obtain the wave crest height and period joint distribution as shown by the black lines in Fig. 3. Please note that in the legend of Fig. 3 the seemingly "grey" line is actually a black line, and it corresponds to the "Level curves enclosing: 99.9". We can obviously notice that the wave crest height and period joint distribution directly obtained from the measured wave elevation data is not symmetrical with respect to the normalized wave period. Therefore, the symmetrical joint distribution model in Longuet-Higgins (1975) is indeed not a suitable model to predict the wave crest height and period joint distributions.
The solid red lines in Fig. 3 represent the empirical joint wave crest height and period distribution as developed by Cavanié et al. (1976). Obviously in this case of Cavanié et al. (1976), the wave crest height and period joint distribution deviates a lot from the corresponding real joint distribution obtained directly from the measured data. [ In Fig. 4, the solid black lines again represent the results of the wave crest height and period joint distribution directly obtained from the measured three-hour wave elevation points at the Poseidon platform. The solid green lines represent the results of the wave crest height and period joint distribution obtained from utilizing our proposed transformed linear simulation method. In applying the transformed linear simulation method, two matrices and as specified in Eq. (21) are first calculated and the elements in these two matrices are directly related to the wave spectrum as shown in Fig. 2 Fig. 2 into a non-Gaussian time history. Next, from this non-Gaussian time history the wave crest height time series are extracted. The time series of the wave crest front period and crest back period are also extracted from this non-Gaussian time history, and adding these two together we can obtain the wave period time series. Then exact Epanechnikov kernel density estimates are carried out to obtain the wave crest height and period joint distribution as shown by the green lines in Fig.  4. We can see that the results from the transformed linear simulation method fit the measured data much better than the results obtained from using the Cavanié et al. (1976) joint distribution model. Moreover, we can obviously notice that the wave crest height and period joint distribution obtained using the transformed linear simulation method is not symmetrical with respect to the normalized wave period. Finally, it is mentioned that in this case it takes only about 16.5 seconds on a desktop computer (HP dx2310MT (VP784PA), Pentium (R) Dual-Core CPU E5300 @2.60 GHz) to obtain the green curves in Fig. 4 by utilizing our proposed transformed linear simulation method.
In Fig. 5, the solid black lines again represent the results of the wave crest height and period joint distribution directly obtained from the measured three-hour wave elevation points at the Poseidon platform. The solid green lines represent the results of the wave crest height and period joint distribution obtained from using the linear simulation method, i.e. by numerically implementing Eqs. (7) and (8) to generate a time series of sea waves with 2000000 elevation points based on the wave spectrum in Fig. 2, and then statistically processing this time series to obtain the wave crest height and period joint distribution. We can obviously find that the results from the linear simulation fit more poorly to the real distribution than the results from the transformed linear simulation do.
In Fig. 6, the solid black lines again represent the results of the wave crest height and period joint distribution obtained directly from the measured three-hour wave elevation points at the Poseidon platform. The solid green lines represent the results of the wave crest height and period joint distribution obtained from using the nonlinear simulation method, i.e. by numerically implementing Eq. (14) to generate a time series of sea waves with 400000 elevation points, and then statistically processing this time series to obtain the wave crest height and period joint distribution. In this case, it takes about 519 seconds on a desktop computer (HP dx2310MT (VP784PA), Pentium (R) Dual-Core CPU E5300 @2.60 GHz) to obtain the green curves in Fig. 6 by using the nonlinear simulation method. If we look more closely on Fig. 6 and Fig. 4 and compare them, we can find that the results from the nonlinear simulation still fit somewhat more poorly to the real distribution than the results from the transformed linear simulation do. For some waves with very high crest heights (e.g. >5 m), the nonlinear simulation method predicts too short corresponding periods. In case of some waves of very long (T c ) and small (A c ) the curves in Fig. 6 do also not agree well. The aforementioned research results demonstrate that our proposed transformed linear simulation method has higher accuracy and efficiency to predict the wave crest height and period joint distribution than the nonlinear simulation method does.

Conclusions
The detailed mathematical formulas of a proposed transformed linear simulation method to predict the wave height and period joint distributions of nonlinear waves have been elucidated in this paper, and the corresponding computer programs have also been developed. The basic points of the methodology proposed are as follows: The proposed transformed linear simulation method is based on a Hermite transformation model in which the transformation is chosen to be a monotonic cubic polynomial, being calibrated such that the first four moments of the transformed model match the moments of the true process. The obtained Hermite transformation model is subsequently used to transform the time history of a linearly simulated wave elevation points into a non-Gaussian time series. Next, from this non-Gaussian time series the wave crest height and wave period time series are extracted. Then exact Epanechnikov kernel density estimates are carried out to obtain the wave crest height and period joint distribution. The proposed new method has been first applied for calculating the wave height and period joint distributions of a sea state with the surface elevation data measured at the Poseidon platform at the coast of Yura, and its accuracy and efficiency have been verified by comparing the results with those obtained from the Cavanié et al. (1976) joint distribution model, from a linear simulation method and from a nonlinear simulation method. The research results in this paper demonstrate that our proposed transformed linear simulation method can be utilized by engineers as a valuable tool to predict the wave height and period joint distributions in their ocean engineering design projects.