Asymptotic Calculation of the Wave Trough Exceedance Probabilities in A Nonlinear Sea

This paper concerns the calculation of the wave trough exceedance probabilities in a nonlinear sea. The calculations have been carried out by incorporating a second order nonlinear wave model into an asymptotic method. This is a new approach for the calculation of the wave trough exceedance probabilities, and, as all of the calculations are performed in the probability domain, avoids the need for long time-domain simulations. The proposed asymptotic method has been applied to calculate the wave trough depth exceedance probabilities of a sea state with the surface elevation data measured at the coast of Yura in the Japan Sea. It is demonstrated that the proposed new method can offer better predictions than the theoretical Rayleigh wave trough depth distribution model. The calculated results by using the proposed new method have been further compared with those obtained by using the Arhan and Plaisted nonlinear distribution model and the Toffoli et al.’s wave trough depth distribution model, and its accuracy has been once again substantiated. The research findings obtained from this study demonstrate that the proposed asymptotic method can be readily utilized in the process of designing various kinds of ocean engineering structures.


Introduction
The wave trough depth is a very important parameter to be considered in the process of designing various kinds of ocean engineering structures. For example, when designing an offshore semi-submersible drilling (or production) platform, it is necessary to define the maximal wave trough depth because the underwater cross-bars of the platform must not be exposed to the air, but at the same time should be sufficiently close to the water surface. The wave trough depth distribution is also very important in the process of designing a tension-leg platform because it is used for the calculation of the tether loads acting on the platform (Toffoli et al., 2008a).
It is known that the wave trough (or crest) distributions obey the theoretical Rayleigh probability law in an ideal Gaussian linear sea (e.g. Romolo and Arena, 2015;Petrova and Soares, 2014;Latheef and Swan, 2013;Longuet-Higgins, 1952). In an ideal Gaussian linear sea model, the individual cosine wave trains linearly superimpose. However, real world waves are generally nonlinear with higher and sharper crests and shallower and smoother troughs (Forristall, 2000). Consequently, the probability distribution of the sea surface elevation in a real ocean tends to deviate from the Gaussian distribution: the larger crest heights are underestimated and the shallower trough depths are overestimated by the Rayleigh distribution. Therefore, the application of the theoretical Rayleigh distribution to the wave troughs (or crests) will become invalid, and other more suitable methods should be applied to predict the distribution of the wave troughs (or crests) for the nonlinear random model of the sea surface elevation.
In the literature, there exist several theoretical and/or empirical models of the wave trough depth distributions of nonlinear random waves. For example, Arhan and Plaisted (1981) and Toffoli et al. (2008a) produced the wave trough probability distribution from the nonlinear Stokes wave model. However, Wang and Xia (2012) have shown that such kinds of theoretical and/or empirical models will sometimes predict the wave characteristic distributions that differ considerably from the true ones. Motivated by the above-mentioned facts, in this paper, a new approach by incorporating a second order nonlinear wave model into an asymptotic method is proposed to calculate the nonlinear trough depth distributions of irregular Stokes waves more accurately and efficiently. The proposed new approach is applied to predict the wave trough exceedance probabilities of a sea state with the surface elevation data measured at the coast of Yura in the Japan Sea. Meanwhile, the wave trough exceedance probabilities of this sea state obtained by using the conventional Rayleigh distribution model and by the Arhan and Plaisted model and Toffoli et al. model are also included for comparison. The calculation results are analyzed, and some findings valuable for engineering design are pointed out.
This paper begins in Section 2 by introducing the knowledge of the irregular Stokes waves. It continues in Section 3 by elucidating the theoretical background of the proposed new approach. In Section 4 the calculation examples and discussions will be provided, with concluding remarks summarized in Section 5.

Irregular Stokes waves
For an ideal linear Gaussian sea, the probability distribution of the wave troughs can be calculated according to the following Rayleigh law: where h t is the height of the wave trough, H s is the significant wave height, and F( ) denotes a probability distribution function. In the Gaussian 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 when the sea states become too severe or as the water depths decrease, the non-linearity of sea waves become more and more relevant, and the ocean surface elevation process deviates significantly from the Gaussian assumption.
For nonlinear Stokes waves, the non-Gaussian effects should be considered when predicting the wave trough distributions. Obviously, the Rayleigh distribution which is good for predicting the trough depths of linear Gaussian waves will become invalid for calculating the wave trough distributions of nonlinear Stokes waves. Arhan and Plaisted (1981) produced a wave trough probability distribution from the nonlinear Stokes wave model, which is given by: where the wave effective steepness R * is given by: where k is the wave number and d is the water depth. Toffoli et al. (2008a) proposed another model for the wave trough exceedance probability as follows: where k is the wave number, and H s is taken to be four times the standard deviation of the measured surface elevation. However, Wang and Xia (2012) showed that such kinds of theoretical and/or empirical models will sometimes predict the wave characteristic distributions that differ considerably from the true ones (For detailed information regarding the differences between the empirical distributions and the true ones, the readers should refer to Section 4.3 and Fig. 15 in Wang and Xia (2012)). Based on the results of extensive direct numerical simulations, Toffoli et al. (2008aToffoli et al. ( , 2008b found that the functions in Eqs. (2) and (4) underestimate the nonlinear wave trough distributions. The deviations are more noticeable at the low probability levels. Furthermore, as the Benjamin-Feir Index (which measures the relative importance of nonlinearity and dispersion) increases, the deviations become even more significant. Motivated by the above-mentioned facts, in this paper, the author proposes a new approach by incorporating a second order nonlinear wave model into an asymptotic method in order to calculate the nonlinear trough depth distributions of irregular Stokes waves more accurately and efficiently.

Proposed new approach
η It has been shown by Wang (2016) that the wave surface elevations (x, t) for the irregular Stokes waves at a specific reference location (say x=0) can be written in the following matrix notation: Please note that a similar expression for the wave surface elevation as Eq. (5) first appeared in Langley (1987). However, Langley (1987 solely dealt with the distribution of the wave surface elevation, which is an entirely different physical quantity from a wave height or a wave trough. Meanwhile the theory behind the proposed new approach in this paper is entirely different from the theory of Langley (1987). In Eq. (5), Q and R are the real symmetric matrices whose nm-th components are and , respectively; s, x and y are the vectors whose n-th components are s n , x n and y n , respectively and . x m and y m are the Gaussian random variables at any fixed time t and have the following properties: By performing an eigenvalue decomposition, Eq. (5) becomes: Λ i where is the diagonal matrix with the eigenvalues in the respective diagonal and P i contains the corresponding eigenvectors per row. Introducing a new set of Gaussian random variables Z j , such that: η we can write the stochastic process as (t) (this is called the Kac-Siegert solution): where Z j denotes 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 .
From Section 2 we know that is the height of the sea level at a fixed point as the function of time t. We here use the so called mean down crossing wave to define a wave, i.e. the wave is considered as a part of a function between the consecutive down crossings of the mean sea level. Assume crossing a mean sea level finite many times and denote the times of down crossings of by , . The trough , say, of the i-th wave is the global minimum of during the interval .
For a specific nonlinear irregular sea state, is a non-Gaussian stochastic process. Here we will use the simplest (but widely effective) model for a non-Gaussian sea where is expressed as the function of a stationary zero mean Gaussian process with variance one ( ), i.e.
Therefore, in order to use the model it is necessary to estimate the transformation G. Here the transformation G relating the non-Gaussian process and the original Gaussian process will be obtained based on the equivalence of the level up-crossing rates of the two processes. Consequently, we can find that in order to calculate the wave trough depth distributions of the non-Gaussian process , the critical task is to calculate the sea level up-crossing rate of . In the following, the principles of a proposed asymptotic method for calculating , will be described.
Let be a strictly non-Gaussian stationary process. For a fixed level u, assuming that is the expected number of times, in the interval [0, 1], the process crosses the u level in the upward direction. We then have the fol-lowing extension of the Rice formula (Machado and Rychlik, 2003): where means that the equality is valid for almost all . In the above equation, z is a realization of and is the derivative of . We see that the computation of requires the estimation of the joint density of and , i. e. .
Although seemingly trivial, the problem of estimating the joint density of and is actually very difficult because for most cases of non-Gaussian processes this density does not have an explicit form. Therefore, in this article we consider a different approach to the approximation of the sea level u up-crossing rate , using the theory of asymptotic expansions. The asymptotic results proposed by Breitung (1988) and further investigated by Hagberg (2004) are employed. In Appendix, we show that the level u up-crossing rate is well approximated by The explicit formula for the constant as well as a justification of this approximation can be found in Appendix. The increasing function is equal to the so-called Hasofer-Linds safety index often employed in the reliability analysis. Note that the estimation of the function is achieved by standard numerical methods. μ(u) η(t) After the level up-crossing rate of the non-Gaussian process is calculated by Eq. (14), it can subsequently be used to obtain the inverse transformation (Wang and Xia, 2013): Next, we turn to the relation between the level up-crossing rate and the probability distribution of the trough depth . As specified in Wang and Xia (2013), we have the following two equations (where h t denotes the nonlinear trough depth): Combining the above two equations, we can obtain the following relationship for the trough depth distribution: , and is calculated by Eq. (15) .

Calculated results and discussions
In this section the proposed new approach is validated based on the information of the measured wave elevation data at the Poseidon platform. The field data considered here were measured in the Japan Sea during the period from 8:12 a. m. on November 24, 1987 to 7:57 a. m. on November 25, 1987. The measurements took place on a floating offshore structure, Poseidon platform, located 3 km off the coast of Yura in the Yamagata prefecture. Fig. 1 is a map showing the position of the coast of Yura. The water depth at the location is approximately 41 m, thus the non-linear effects (finite water depth effects) are expected to affect the surface elevation. The surface elevation data were measured with three ultrasonic wave gauges located on the sea floor and three sets of data each consisting of 85547 wave elevation points had been obtained during the measurement. The sampling frequency was 1 Hz.
Because a sea state remains stationary only within a short time period of less than three hours (see, e.g. page 197 of Chakrabarti (2005)), we should divide each of the above three sets of data into eight short records respectively when we calculate the wave trough distributions. In this study we have chosen a three-hour record of wave elevation points time series measured with the #1 ultrasonic wave gauge (Data source is from http://www.maths.lth.e/matstat/wafo/ download/index.html. Please note that these data are used for academic research only and not for commercial purposes). The significant wave height of this three-hour record of wave elevation points time series is H s = 5.21 m, and the spectral peak period is T p = 13.13 s. Fig. 2 shows the wave spectrum corresponding to this specific three-hour record of wave elevation points time series. Fig. 3 shows the first 450 wave elevation points in this three-hour time series measured with the #1 ultrasonic wave gauge. We can obviously find that the waves are nonlinear. Fig. 4 shows the calculated wave trough exceedance probabilities for this measured sea state. The solid green line in Fig. 4 represents the results of the wave trough ex-ceedance probabilities directly obtained from the measured Poseidon data (10700 wave elevation points). For obtaining these results, the wave trough time series were first extracted from these 10700 wave elevation points. Then exact Epanechnikov kernel density estimates were carried out for obtaining the probability density function of the heights of the wave troughs. Next, the cumulative trapezoidal numerical integration was performed on the above-mentioned probability density function for getting the probability distribution (F) of the depths of the wave troughs. Finally, the wave trough exceedance probabilities were obtained by using the formula P=1-F. These calculation results based on the measured data are used as the benchmark against which the accuracy of the results from the asymptotic method and    from the existing trough distribution models is checked. In Fig. 4, the solid blue line represents the results of the wave trough exceedance probabilities obtained from utilizing the proposed asymptotic method. In applying the asymptotic method, the G -1 function is calculated by using Eqs. (14) and (15) in Section 3 and the equations in the Appendix. We can see that the results from the asymptotic method fit the results from the measured data fairly well.
In Fig. 4, the small blue "•" represents a wave trough exceedance probability predicted by the Arhan and Plaisted model in which the value of H s was taken to be four times the standard deviation of the measured surface elevation. We can find that the Arhan and Plaisted model gave poorer predictions in the wave trough region of about [2.8 m, 4.3 m] than the asymptotic method did. In Fig. 4, the solid red line represents the wave trough exceedance probabilities calculated by using a theoretical Rayleigh model. We notice that the theoretical Rayleigh method overestimates the exceedance probabilities in the wave trough region of about [1.6 m, 4.3 m].
In Fig. 4, the solid pink line represents the wave trough exceedance probabilities predicted by using the Toffoli et al. (2008a) model, i.e. by using the following equation in which the value of H s was taken to be four times the standard deviation of the measured surface elevation.
where k is the wave number. We can find that the Toffoli et al. (2008a) model gave poorer predictions in the wave trough region of about [2.9 m, 4.3 m] than the asymptotic method did. If we look at Fig. 4 more closely, we can find that it is true that the presented asymptotic wave trough distribution is the closest to the result from the measured Yura data. But it still slightly under predicts the results. Here some explanations for the difference are given as follows.
The asymptotic analysis method in the present paper is a numerical method. The functional transformation in the present paper is related to the level up-crossing rates of a non-Gaussian process, and these level up-crossing rates are calculated by using an asymptotic analysis method. From Appendix we can find that the essence of the asymptotic analysis method in the present paper is to use the second order reliability analysis method to calculate the up-crossing rates for a noncentral Chi square process. Owing to the curse of dimensionality in the probability-of-failure calculation, numerous methods are used to simplify the numerical treatment of the integration process. The Taylor series expansion is often used to linearize the limit-state g(U)=0. In the present approach, the second-order Taylor series expansion is used to perform the second-order approximation of the response surface g(U) =0 at the most probable failure point. This approach is actually a mathematical optimiza-tion problem for finding the point on the response surface that has the shortest distance from the origin to the surface in the standard normal space.
As can be noticed from the aforementioned explanations, the proposed asymptotic analysis method still involves an approximation procedure, and this can lead to the difference between the predicted asymptotic distribution and the measured Yura data.
Finally, it is mentioned that it took less than 2 seconds on a desktop computer to obtain the results of the wave trough exceedance probabilities represented by the blue solid line in Fig. 8 by applying the asymptotic method. The accuracy and efficiency of the asymptotic method for calculating the wave trough distributions of irregular Stokes waves can thus be validated.

Concluding remarks
The detailed mathematical procedures and formulas of an asymptotic method for calculating the wave trough distributions of irregular Stokes waves have been elucidated in this article, and the second order Stokes wave model has been integrated into this asymptotic method. The proposed asymptotic method has been applied for calculating the wave trough exceedance probabilities of a sea state with the surface elevation data measured at the Poseidon platform. It is demonstrated that the asymptotic method can offer better predictions than those from using the theoretical and/or empirical wave trough distribution models. The research findings obtained from this study demonstrate the suitability of utilizing the proposed asymptotic method for engineering purposes.  (4), 3022-3042. Toffoli, A., Bitner-Gregersen, E., Onorato, M. and Babanin, A.V., 2008a. Wave crest and trough distributions in a broad-banded directional wave field, Ocean Engineering, 35(17-18), 1784-1792. Toffoli, A., Onorato, M., Bitner-Gregersen, E., Osborne, A.R. and Babanin, A.V., 2008b

Appendix:
In this paper we are interested in the measurements of the sea surface elevation at a fixed origin point, hence we may write η (t) = X (t) . (A1) The following derivation process follows that in Baxevani et al. (2005Baxevani et al. ( , 2009): With comparison with Eq. (23) we may write the process in the form (Baxevani et al., 2005(Baxevani et al., , 2009): where is a vector-valued stationary Gaussian process, such that for each t, and the variables and are independent (j and k are two different integers, , ). Let us also denote by , the derivative of vector . Then the joint density of the vectors and is normal with , where (Baxevani et al., 2005(Baxevani et al., , 2009 while I is the identity matrix (note that the matrices , need not be identity matrices). Note that Z 0 (t) = X 0 (t)/ β 0 , where is the Gaussian part of the process and , is independent of the processes Z j (t), j=1, 2,..., N. (Obviously and ). The linear part is a Gaussian process with a spectrum S η (ω) and variance in which the equality is for going to infinity. The quadratic correction term X q (t) = has mean zero and variance obtained as . By the independence for fixed t of the different and because and are uncorrelated, the variance of is the sum of the variances of the terms in Eq. (A2) (Baxevani et al., 2005(Baxevani et al., , 2009 The following generalization of Breitung's approximation (Breitung, 1988) can be found in Baxevani et al. (2005Baxevani et al. ( , 2009 Let g: be a function such that the surface has a point such that and for all other . By we denote both the vector and the column matrix. Suppose is an n-dimensional, stationary, differentiable, Gaussian vector process, and let denote its derivative. The correlation of the vector is denoted by , , β < 0 For a family of processes , under some mild technical assumptions, the intensity of zero upcrossings is given by (for detailed derivation process, please refer to Belyaev (1968), Breitung (1988), andHagberg (2004)): In order to perform calculations for the general case of a second order sea, we have to construct somewhat artificial asymptotics. The idea is as follows. Fix the level u, assumed to be large, and let: where b is a column vector containing and Λ = diag[γ 1 , γ 2 , γ 3 , ......, γ n ] is the diagonal matrix, with . Assume that there is only one point on the surface β u := ∥x u ∥ of minimal distance to the origin, and define . Let The process crosses the level 0 when the process crosses the level u. Hence, with , for the specific level u, we have equal to the u-upcrossing rates for the process . Therefore, it is reasonable to believe that if is large, then the term for is small. Hence the approximation of omitting the term in Eq. (A6) is suitable.
β u Note that, here we indicate c's dependence of . Evaluating the terms gives Now we can use the following approximation: It should be noted that the point of minimum norm, can be found by standard optimization methods. After finishing the critical task of calculating the level up-crossing rate using Eq. (A12), we can then calculate the wave crest height distributions of the non-Gaussian process .