Estimation of Wave Crest Amplitudes Distribution and Freak Wave Occurrence in A Short Crested Mixed Sea

In this study we have for the first time proposed a novel transformed linear simulation method for the estimation of wave crest amplitudes distribution and freak wave occurrence in a short crested mixed sea with a bimodal 3D spectrum. For implementing the proposed transformed linear simulation method, a Hermite transformation model expressed in a monotonic cubic polynomial has been constructed so that the first four moments of the original true process match the corresponding moments of the transformed model. The proposed novel simulation method has been applied to forecast the freak wave occurrence in two short crested mixed sea states, one with a directional wave spectrum based on the measured surface elevation data at the coast of Yura, and the other one with a typical directional bimodal Torsethaugen wave spectrum. It is shown in the two cases that the proposed novel simulation method can offer more accurate forecasting results than those obtained from the traditional linear simulation method or by using Rayleigh distribution model. It is also demonstrated in this article that the proposed novel simulation method is more efficient than the nonlinear simulation method.


Introduction
Freak waves are events that feature a single, unusually large and steep crest, and post severe hazards to ships and ocean engineering structures. Freak waves are not necessarily the biggest waves found on the water; they are, rather, unusually large waves for a given sea state. Once considered mythical and lacking hard evidence for their existence, freak waves are now proven to exist and known to be a natural ocean phenomenon. Eyewitness accounts from mariners and damage inflicted on ships have long suggested their occurrence. The first scientific evidence of freak waves came with the recording of a freak wave by the Gorm platform in the central North Sea in 1984. A stand-out wave was detected with a wave crest height of 11 m in a relatively low sea state (Sand et al., 1990). However, the wave that caught the attention of the scientific community was the digital measurement of the "Draupner wave", a freak wave at the Draupner E platform in the North Sea on January 1, 1995, with a maximum wave height of 25.6 m (crest peak elevation of 18.5 m). During that event, minor damage was also inflicted on the platform, far above sea level, confirming that the reading was valid (Haver, 2003). Hayer and Andersen (2000) collated evidence that freak waves were not the rare realizations of a typical or slightly non-Gaussian sea surface population but rather they were the typical realizations of a rare and strongly non-Gaussian sea surface population of waves. An interesting explanation was given earlier by Dean (1990) to indicate that both nonlinearity and directionality are primary possible causes of freak waves. Laboratory studies have demonstrated that freak waves can be generated through non-linear wave-wave interactions in a two-dimensional wave flume (Stansberg, 1990). Mori et al. (2002) showed that a well defined freak wave may occur in the developed wind-wave condition with single peak directional spectra. The crest and trough amplitude distributions of the observed sea waves including freak waves are different from the Rayleigh distribution. Mori et al. (2002) further concluded that the exceedance probabilities of the crest amplitude of the freak wave are underestimated by the Rayleigh distribution. Therefore, accurate calculation of the exceedance probabilities of the wave crest amplitudes is critical for the prediction of the occurrence of freak waves in a nonlinear (or non-Gaussian) sea.
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. Based on this specific wave spectrum, there are several approaches to calculate the exceedance probabilities of the wave crest amplitudes. The first and most straightforward way is to resort to an empirical (or theoretical) model. Chakrabarti (2005) mentions the Rayleigh model of wave crest amplitude distribution of random ocean waves. However, the previous research work of Wang and Xia (2013) has shown that the Rayleigh model systematically underestimates the wave crest amplitude distribution of nonlinear random ocean waves. From a specific wave spectrum, there is another way to calculate the exceedance probabilities of the wave crest amplitudes. Based on the wave spectrum at the Gullfaks C platform, Wang and Xia (2012) calculated the wave crest amplitude distribution using a linear simulation method. However, their calculation results showed that the linear simulation method will predict overly non-conservative probability distributions of the wave crest amplitudes in an actually nonlinear sea state. This will result in the design of unsafe ocean engineering structures. Wang and Xia (2012) used a nonlinear simulation method to compute the wave crest amplitude distribution for the nonlinear irregular waves. They had verified that their nonlinearly simulated wave crest amplitude distribution is more accurate than that calculated by using the linear simulation method. 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 ocean engineering projects. Wang (2014) proposed a Transformed Rayleigh method for calculating the wave crest amplitude distribution of shallow water nonlinear waves. It was demonstrated by Wang (2014) that the Transformed Rayleigh method has a higher efficiency and equivalent accuracy, in comparison with the nonlinear simulation method. However, the Transformed Rayleigh method developed by Wang (2014) is only applicable to a sea state with an ideal 2D (long crested) wave spectrum, and this drawback will certainly hinder its usefulness in solving some real world ocean engineering problems. The prediction method of the wave crest amplitude distribution presented by Wang (2014) assumes that wave energy is traveling in a specific direction, commonly considered the same direction as the wind. In this respect, the wave spectrum in Wang (2014) may be considered as a long crested spectrum or a uni-directional spectrum. In reality, however, wind-generated wave energy does not necessarily propagate in the same direction as the wind; instead, the energy usually spreads over vari-ous directions. Thus, for an accurate description of random seas, it is necessary to clarify the spreading status of energy. The wave spectrum representing energy spreading over multi-directions is called a multi-directional spectrum or a short crested (3D) spectrum. Information on wave directionality is extremely significant for the design of marine systems such as ships and ocean engineering structures. This is not only because the responses of a system in a seaway computed using a unidirectional (2D long crested) wave spectrum are overestimated; the associated coupled responses induced by waves from other directions are also disregarded.
Furthermore, all the above-mentioned wave spectra (no matter a 2D spectrum or a 3D spectrum) are unimodal wave spectra. It is well known that not all sea states have unimodal wave spectra and narrow (or finite) spectral bandwidth. Frequently, sea states are due to the coexistence of various wave systems. In particular, local wind waves often develop in the presence of some background low frequency swell coming from distant storms, and the resulting mixed sea states will have bimodal wave spectra (Wang and Xia, 2013).
Therefore, it is of vital practical importance to develop an accurate and efficient method for predicting the wave crest amplitudes distribution and freak wave occurrence in a short crested mixed sea.
Motivated by the afore-mentioned facts, the present study will first propose a novel transformed linear simulation method for the estimation of wave crest amplitudes distribution and freak wave occurrence in a short crested mixed sea with a bimodal 3D spectrum. For implementing the proposed transformed linear simulation method, a Hermite transformation model expressed in a monotonic cubic polynomial will be constructed so that the first four moments of the original true process match the corresponding moments of the transformed model. This novel simulation method will be utilized for calculating wave crest amplitude exceedance probabilities of two sea states, one with a directional (3D) wave spectrum based on the measured surface elevation data at the coast of Yura, and the other one with a typical directional (3D) Torsethaugen wave spectrum. The calculated exceedance probabilities of a specific critical wave crest amplitude are directly related to the occurrence of freak waves. It will be demonstrated that the novel simulation method can produce predictions more accurate than those obtained from traditional linear simulation method and from Rayleigh crest amplitude distribution model. It will also be shown that the efficiency of the novel simulation method is higher than that of the nonlinear simulation method.

Nonlinear random waves and the transformed linear simulation method
Waves in an idealized linear Gaussian random sea have crest-trough symmetry. However, it is known that real η (x, t) world ocean surface elevation process deviates from the Gaussian assumption, i.e. the wave crests are becoming steeper and higher and the wave troughs are becoming flatter and shallower than expected under the Gaussian assumption. If the longitudinal coordinate is denoted by x and time is denoted by t, the free surface elevation for a nonlinear random sea state can then be written as: in which are the first-order linear components and are the second-order nonlinear correction compon- (1), usually is chosen to be a sufficiently large positive integer, denotes a random complex valued wave amplitude that can be calculated based on a specific wave spectrum and the angular frequency . For a short-crested sea state with a multi-directional wave spectrum in which denotes the wave direction, an equivalent frequency spectrum can be obtained by taking the directional (3D) wave spectrum and integrating the energy over all directions to give the total energy at each frequency. Furthermore, in Eq. (1) is a specific wave number that is related to the through the dispersion relation, and is the uniformly distributed phase angle. Finally, the terms and in Eq.
(1) are secondorder transfer functions that can be calculated (for a constant water depth d) by using the following equations: ; (2) .
(1) will lead to linearly simulated wave time series, and this is called the linear simulation method. Implementing Eqs.
(1)-(3) altogether will lead to nonlinearly simulated wave time series, and this is called the nonlinear simulation method. However, because Eqs. (2) and (3) contain N 2 sum frequencies components and another N 2 difference frequencies components, the nonlinear simulation based on Eqs. (1)-(3) will become very time consuming when N is large. In order to improve the simulation efficiency for generating shallow water nonlinear waves, in this article we propose a novel transformed linear simulation method whose theoretical background is explained as follows. η Without loss of generality, we set x=0 in and simplify it as . The non-Gaussian wave elevation process can be modeled as a function of a standard Gaussian process : In Eq. (4), G (·) is a continuously differentiable function with positive derivative. The transformation G (·) performs the appropriate nonlinear translation and scaling so that is always normalized to have mean zero and variance one. In the present study, a Hermite transformation model expressed in a monotonic cubic polynomial will be constructed for G (·) so that the first four moments of the original true process match the corresponding moments of the transformed model, i. e.
In Eq. (5), and are the mean and standard deviation values of the second-order nonlinear random waves. For N=4 moments, the expressions for the coefficients and are (see Wang (2014)): The coefficients and in Eqs. (6) and (7) are the skewness and kurtosis of the surface elevation process respectively. From Eqs. (5)-(7) we can see that for calculating G(·) function, we need to know the values of , , , and of the second-order nonlinear random waves. For a chosen sea state with a specific wave spectrum , we can obtain the values of the coefficients and according to the theories of Langley (1987). Then the values of , , , and can be calculated by using the following equations (Langley, 1987):

Calculation examples and discussions of the calculation results
In the following we will validate the accuracy and efficiency of the proposed transformed linear simulation method by utilizing several calculation examples. Specifically, the proposed transformed linear simulation method will be applied for calculating the wave crest amplitude exceedance probabilities of two specific sea states, one with a directional (3D) wave spectrum based on the measured surface elevation data at the coast of Yura, and the other one with a typical directional (3D) Torsethaugen wave spectrum. The calculated exceedance probabilities of specific critical wave crest amplitude are directly related to the occurrence of freak waves. The measured surface elevation data at the coast of Yura were obtained at the location 3 km off Yura fishing harbor facing the Sea of Japan. The observations were carried out during the period from 2:11 a. m. to 5:11 a. m. on November 25, 1987 by the Ship Research Institute, Ministry of Transport of Japan. Temporal sea surface elevations were measured with ultrasonic-type wave gages installed at three points in 42 m water depth. The sampling time interval during the measurement was 1s. Based on these measured Yura coast surface elevation data, a directional (3D) wave spectrum was estimated by using the Maximum Likelihood Method and is shown in Fig. 1. We can obviously notice that this wave spectrum is bimodal, indicating that the corresponding sea state is a short crested mixed sea state. Fig. 2 shows the measured wave elevation time series by the wave gage at the origin (x=0, y=0) and this time series contains 10800 wave elevation points. is larger than 7.8246 m (twice the significant wave height, i.e., 2H s = 7.8246 m). We have carefully checked and have found that in the time series in Fig. 2 there are no other waves whose heights are twice greater than the significant wave height (H s ), i.e., there are no other freak waves in the time series in Fig. 2. At first glance it seems that there is another freak wave at the time instant of about 5700 s, however, the zoom-in plot in Fig. 4 reveals that this wave has a too small wave trough depth, and therefore the total height of this wave is not great although it has a very large   crest height, i.e. this wave is not a freak wave.
In a real world ocean engineering project, the known information regarding a sea state is typically a specific wave spectrum rather than a measured wave elevation time series. We then need to forecast the occurrence probability of the extreme waves (including freak waves) in this sea state based only on the information of the wave spectrum so that we can design a safe ocean engineering structure. In the Introduction part of this paper, we have mentioned that based on a specific wave spectrum, an empirical (or theoretical) model or numerical simulation methods can be used to calculate the exceedance probabilities of the wave crest amplitudes (which is directly related to the occurrence probability of the extreme waves (including freak waves)). In the following, taking the directional (3D) wave spectrum in Fig. 1 as a calculation example, we will propose a novel transformed linear simulation method to calculate the exceedance probabilities of the wave crest amplitudes, and compare its accuracy and efficiency with those of a theoretical model or the linear and nonlinear simulation methods. Fig. 5 shows our calculated wave crest amplitudes exceedance probabilities based on the directional (3D) wave spectrum in Fig. 1 and the measured Yura coast wave data. In Fig. 5 the black dots show the calculation results of the wave crest amplitudes exceedance probabilities directly obtained from the measured Yura coast wave data that contains 10800 wave elevation points. These black dots based on the measured Yura coast wave data are used as the benchmark against which the accuracy of the results from various numerical simulation methods and from an existing wave crest amplitudes model is checked. In Fig. 5 the solid red line shows the results of the wave crest amplitudes exceedance probabilities obtained from using the Rayleigh distribution model as expressed in the following equation: A c h H s In Eq. (12) represents the wave crest amplitude, is the specific value of the wave crest amplitude, is the significant wave height. We can clearly see that the wave crest amplitudes exceedance probabilities obtained by using the Rayleigh distribution model deviate a lot from the corresponding benchmark results directly obtained from the measured Yura coast wave data. This is not a surprise because the Rayleigh distribution is an ideal linear model and the measured Yura coast wave elevation series are obviously nonlinear. In Fig. 5 the solid green line represents the results of the wave crest amplitude exceedance probabilities obtained by using the linear simulation method based on the directional (3D) wave spectrum in Fig. 1. The linear simulation process started with taking the directional (3D) wave spectrum in Fig. 1 and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a wave elevation time series of 8.0×10 6 points by applying Eq. (1) omitting the term . Next the wave crest amplitudes time series were extracted from the 8.0×10 6 wave elevation points. Then exact Epanechnikov kernel density estimates were carried out for obtaining the probability density function of the amplitudes of wave crests. Next, cumulative trapezoidal numerical integration was performed on the above-mentioned probability density function for getting the probability distribution (F) of the amplitudes of wave crests. Finally, the wave crests amplitudes exceedance probabilities were obtained by using the formula P=1F. It took about 9.5 seconds to finish the entire linear simulation process as mentioned above on a desktop computer (Lenovo ThinkCenter, Intel (R) Core(TM) i7-4790 CPU@3.60 GHz). We can find that the wave crest amplitudes exceedance probabilities obtained by using the linear simulation method also deviate a lot from the corresponding benchmark results directly obtained from the measured Yura coast wave data.
From Fig. 5 we can read (actually from a detailed calculation report) the wave crest amplitude (h=4.108 m) exceedance probabilities obtained using various methods and these values are summarized in Table 1. From Table 1 we   5. Comparison between the wave crest amplitudes exceedance probabilities from the linear simulation, from the nonlinear simulation, from the transformed linear simulation, from the Rayleigh distribution and from the wave crest amplitudes exceedance probabilities results from the measured Yura coast wave data.
can find that the wave crest amplitude (h=4.108 m) exceedance probability (P=0.0001) obtained by using the linear simulation method is much smaller than the corresponding benchmark exceedance probability (P=0.0013) directly obtained from the measured Yura coast wave data. A 4.108m wave crest amplitude h=4.108 m) corresponds to a wave height (H) of 7.8246 m. The reason is that the wave height is equal to the sum of the wave crest amplitude and the wave trough depth. Meanwhile, Wang and Wang (2018) has shown that in a nonlinear sea a wave crest amplitude (h) is not equal to half of the wave height (H/2), but rather is approximately equal to 1.05×(H/2). This is why the 4.108 m wave crest amplitude (h=4.108 m) corresponds to the wave height (H) of 7.8246 m. From the previous paragraphs we know that the significant wave height of the measured wave elevation time series in Fig. 2 is H s = 3.9123 m. Therefore 7.8246 m is twice the significant wave height of the measured wave elevation time series at the coast of Yura shown in Fig. 2. (i.e. 2H s =7.8246 m). Therefore, the wave crest amplitude (when h=4.108 m) exceedance probability corresponds to the exceedance probability of the wave height which is twice the significant wave height. Therefore, in this case the wave crest amplitude (when h=4.108 m) exceedance probability is a critical value indicating that freak waves have occurred in a short crested mixed sea. Consequently, it is of uttermost importance to accurately calculate this wave crest amplitude (when h=4.108 m) exceedance probability so that we can design safer ocean engineering structures to withstand the possible impacts of freak waves.
In order to more accurately calculate the wave crest amplitudes exceedance probabilities we had tried to use a nonlinear simulation method. In Fig. 5 the pink "+" represents the result of the wave crest amplitude exceedance probability obtained using the nonlinear simulation method based on the directional (3D) wave spectrum in Fig. 1. The nonlinear simulation process started with taking the directional (3D) wave spectrum in Fig. 1 and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a nonlinear wave elevation time series of 800000 points by applying Eqs. (1)-(3). Further mathematical and statistical processing finally leads to the wave crests amplitudes exceedance probabilities by the non-linear simulation method as shown in Fig. 5. It took about 314 seconds to finish the entire nonlinear simulation process as mentioned above on a desktop computer (Lenovo ThinkCenter, Intel(R)Core(TM)i7-4790 CPU@3.60 GHz). We can find that in the wave crest amplitudes region [0 4.108 m] the exceedance probabilities obtained by using the nonlinear simulation method fit quite well with the corresponding benchmark results directly obtained from the measured Yura coast wave data.
From Fig. 5 we can once again read the wave crest amplitude (h=4.108 m) exceedance probability obtained by using the nonlinear simulation method and this value is also summarized in Table 1. From Table 1 we can find that the wave crest amplitude h=4.108 m) exceedance probability (P=0.00063) obtained by using the nonlinear simulation method does not have much difference from the corresponding benchmark exceedance probability (P=0.0013) directly obtained from the measured Yura coast wave data. As explained before, in this case the wave crest amplitude (when h=4.108 m) exceedance probability is a critical value indicating that freak waves have occurred in a short crested mixed sea. Therefore, the more accurate wave crest amplitude (when h=4.108 m) exceedance probability predicted by utilizing the nonlinear simulation method will obviously help us to design safer ocean engineering structures. However, we have noticed that in this case it took about 314 s to finish the entire nonlinear simulation process. Although this efficiency is practically acceptable, we can still pursue other more efficient and accurate methods.
In order to further improve the simulation efficiency, we had tried to use a transformed linear simulation method to more efficiently and accurately calculate the wave crest amplitudes exceedance probabilities. In Fig. 5 the solid blue line represents the results of the wave crest amplitudes exceedance probabilities obtained by using the transformed linear simulation method based on the directional (3D) wave spectrum in Fig. 1. The transformed linear simulation process also started with taking the directional (3D) wave spectrum in Fig. 1 and integrating the energy over all directions to give the total energy at each frequency. The obtained equivalent frequency spectrum was then utilized to generate a linear wave elevation time series of 8000000 points. Then the theories in Eqs. (4)-(11) were utilized to calculate a transformation function (G(·) function) and this G(·) func- tion was then used to transform the aforementioned linearly simulated wave elevation time series of 8000000 points into an equivalent "nonlinear" time series. Further mathematical and statistical processing finally leads to the wave crests amplitudes exceedance probabilities as represented by the solid blue line. It took about 10.5 seconds to finish the entire transformed linear simulation process as mentioned above on a desktop computer (Lenovo ThinkCenter, Intel(R)Core(TM)i7-4790CPU@3.60 GHz). We can find that in the wave crest amplitudes region [0 4.108 m] the exceedance probabilities obtained by using the transformed linear simulation method fit quite well with the corresponding benchmark results directly obtained from the measured Yura coast wave data. From Fig. 5 we can once again read the wave crest amplitude (h=4.108 m) exceedance probability obtained using the transformed linear simulation method and this value is also summarized in Table 1. From Table 1 we can find that the wave crest amplitude (h=4.108 m) exceedance probability (P=0.000642) obtained using the transformed linear simulation method does not have much difference from the corresponding benchmark exceedance probability (P=0.0013) directly obtained from the measured Yura coast wave data, and this predicted probability value (P=0.000642) is slightly more accurate than the exceedance probability value (P=0.00063) obtained using the nonlinear simulation method. As explained before, in this case the wave crest amplitude (when h=4.108 m) exceedance probability is a critical value indicating that freak waves have occurred in a short crested mixed sea. Therefore, the more accurate wave crest amplitude (when h=4.108 m) exceedance probability predicted by utilizing the transformed linear simulation method will surely help us to design safer ocean engineering structures in a more efficient way.
In our research we have also validated the accuracy and efficiency of our proposed transformed linear simulation method for predicting the wave crest amplitudes distribution and freak wave occurrence by using another set of measured Yura coast wave data. However, in order to shorten our paper, these validation results are not included in this paper any more.
We have demonstrated the accuracy and efficiency of the transformed linear simulation method in predicting the crest amplitudes distribution and freak wave occurrence in a special case when both a measured wave data record and a directional (3D) wave spectrum are available for a sea state. However, in most other ocean engineering projects that we will meet in the real world, the known information regarding a sea state usually will only contain a specific wave spectrum. That is to say, we usually will not have a measured wave data record for the sea state in these engineering projects due to economic or other reasons. In these general cases, the transformed linear simulation method also undoubtly has higher accuracy and efficiency in predicting the crest amplitudes distribution and freak wave occurrence, and in the following we will carry out a specific calculation example to demonstrate this fact.
In our specific example, we will calculate the wave crest amplitude exceedance probabilities for a specific sea state with a 3D bimodal Torsethaugen spectrum (as shown in Fig. 6) with a significant wave height =5.2 m, a spectral peak period =7 s, a peakedness factor =0.5, and a cosine squared spreading function with the spreading parameter equal to 15. The water depth for this short crested mixed sea state is taken to be 42 m. Fig. 7 shows our calculated wave crest amplitudes exceedance probabilities based on the directional (3D) Torsethaugen wave spectrum in Fig. 6. In Fig. 7 the solid blue line shows the calculation results of the wave crest amplitudes exceedance probabilities obtained by utilizing the transformed linear simulation method (in this process a wave elevation time series of 8000000 points were generated). It took about 11 seconds to finish the entire transformed linear simulation process as mentioned above on the same desktop computer as before). In Fig. 7 Fig. 7. Comparison between the wave crest amplitude exceedance probabilities from the linear simulation based on the 3D Torsethaugen spectrum, from the nonlinear simulation based on the 3D Torsethaugen spectrum, from the Rayleigh distribution and from the wave crest amplitude exceedance probabilities results from the transformed linear simulation based on the 3D Torsethaugen spectrum. ceedance probabilities obtained by using the Rayleigh distribution model as expressed in Eq. (12).
In Fig. 7 the solid green line represents the results of the wave crest amplitude exceedance probabilities obtained by using the linear simulation method (in this process a wave elevation time series of 8.0×10 6 points were generated). In Fig. 7 a pink + represents the result of a wave crest amplitude exceedance probability obtained by using the nonlinear simulation method (in this process a wave elevation time series of 8.0×10 6 points were generated). It took about 299 seconds to finish the entire nonlinear simulation process to obtain the wave crest amplitudes exceedance probabilities on the same desktop computer as before. We can find that the wave crest amplitudes exceedance probabilities obtained by using the transformed linear simulation method fit quite well with the corresponding probabilities values obtained by using the nonlinear simulation method. This fact can substantiate that using the transformed linear simulation method we can efficiently obtain wave crest amplitudes exceedance probabilities values as accurate as those obtained by using the time-consuming nonlinear simulation method.
From Fig. 7 we can read the wave crest amplitude (h=5.46 m) exceedance probabilities obtained using the various methods and these values are summarized in Table 2. From Table 2 we can find that the wave crest amplitude (h=5.46 m) exceedance probability (P=0.00086) obtained by transformed linear simulation method has only a 1.2% difference in comparison with the exceedance probability (P=0.00085) value obtained via the nonlinear simulation method. A 5.46-m wave crest amplitude (h=5.46 m) corresponds to a wave height (H) of 10.4 m. From the previous paragraphs we know that 10.4 m is twice the significant wave height for the directional (3D) Torsethaugen wave spectrum in Fig. 6. (i.e. 2H s =2×5.2 m=10.4 m). Therefore, in this calculation example the wave crest amplitude (when h=5.46 m) exceedance probability corresponds to the exceedance probability of a wave height which is twice the significant wave height. Therefore, in this case the wave crest amplitude (when h=5.46 m) exceedance probability is a critical value indicating that freak waves have occurred in a short crested mixed sea state with a directional (3D) Torsethaugen spectrum in Fig. 6. Consequently, it is of uttermost importance to accurately calculate this wave crest amplitude (when h=5.46 m) exceedance probability so that we can design safer ocean engineering structures. Therefore, in this calculation example, we have demonstrated the accuracy and efficiency of the transformed linear simulation method in predicting the crest amplitudes distribution and freak wave occurrence in a general case when only a directional (3D) wave spectrum is available for a short crested mixed sea state.

Conclusions
In this study we have for the first time proposed a novel transformed linear simulation method for the estimation of wave crest amplitude distribution and freak wave occurrence in a short crested mixed sea with a bimodal 3D spectrum. For implementing the proposed transformed linear simulation method, a Hermite transformation model expressed in a monotonic cubic polynomial has been constructed so that the first four moments of the original true process match the corresponding moments of the transformed model. Firstly, in our study we have demonstrated the accuracy and efficiency of the proposed transformed linear simulation method in predicting the crest amplitudes distribution and freak wave occurrence in a special case when both a measured wave data record and a directional (3D) wave spectrum are available for a sea state. In this special case we have demonstrated the accuracy and efficiency of the proposed transformed linear simulation method in calculating the wave crest amplitudes exceedance probabilities. The calculated exceedance probabilities of a specific critical wave crest amplitude are directly related to the occurrence of freak waves. It is demonstrated in this special case that the proposed transformed linear simulation method can offer more accurate predictions than those by using the traditional linear simulation method and by using the Rayleigh crest amplitude distribution model. It is also shown in this special case that the efficiency of the proposed transformed linear simulation method is higher than that of the nonlinear simulation method.
However, for most other practical projects that people will meet in the real world, the known information regarding a sea state usually will only contain a specific wave spectrum. That is to say, people usually will not have a measured wave data record for the sea state in these engin- eering projects due to economic or other reasons. In these general cases the proposed transformed linear simulation method will also have higher accuracy and efficiency in predicting the crest amplitudes distribution and freak wave occurrence, and in our paper we have carried out a specific calculation example to demonstrate this fact. In this specific calculation example, we have demonstrated the accuracy and efficiency of the proposed transformed linear simulation method in calculating the wave crest amplitudes exceedance probabilities of a sea state with a directional (3D) bimodal Torsethaugen wave spectrum. The research results in this paper demonstrate that the proposed transformed linear simulation method can help people to more accurately and efficiently predict the freak wave occurrence in a short crested mixed sea so that they can design safer ocean engineering structures. Finally, from Fig. 5 we can notice that there is still a need for improving the developed method so that we can even more accurately calculate the exceedance probabilities of very rare, exceptionally high waves. In the future we will pursue further research along this direction and establish prediction methods that also take into account the effect of third-order nonlinearities of the sea waves.