Comparative Study on Numerical Calculation of 2-D Random Sea Surface Based on Fractal Method and Monte Carlo Method

Based on fifty one groups of data on direction distribution measured from buoy sites, several important spectrum parameters including distribution characteristics of the measured data’s spectrum, the Wen’s direction spectrum and the Donelan function are analyzed from the perspectives of standard deviation of directional distribution function and statistical results. Then, a numeric calculation model based on the Monte Carlo method is proposed in this work. At the same time—based on Weierstrass self-affine fractal function—numeric simulation of random sea surface is conducted by modifying the bilateral power law. The analysis of the numeric calculation results under different wind directions, speeds and fetches proves that both methods can be adopted for different water directional distributions and target spectrum models. In addition, this study compares the characteristic wave within different distribution frequency domains in the main wave direction and in its orthogonal direction. It is proved that the fractal method cannot fully reflect the anisotropy of gravity wave and tension wave in the superposition direction, however, it can maintain the similarity of overall energy part with the rough length of the waves. Moreover, there are still limitations of the fractal method in the numeric modeling of undeveloped sea surface.


Introduction
In the physical ocean domain, the sea spectrum is defined as the function distribution spectrum of power density on the sea surface. The sea spectrum is calculated by means of the Fourier transform using the correlation function of height fluctuation on the sea surface, which can reflect the statistical distribution characteristics of the energy contained in the wave in terms of propagation direction and wavelength [1]. Thus, the sea spectrum is also one of the most basic methods used to describe the sea surface. In the field of ocean wave numeric calculation, the sea surface is usually regarded as a superposition of a large number of sine and cosine waves, which usually have different propagation directions and wavelengths [2][3][4]. The spectrum model describes the distribution characteristics of the power spectral density of ocean waves in terms of different propagation directions and wavelength ranges. According to the growth state of the sea, waves can be divided into two types: steady-state spectral model and unsteady-state spectral model. In fluid mechanics, gravity wave refers to a wave in a liquid medium or the interface between two media and its restoring force comes from gravity or buoyancy. During the wave motion, the restoring force required for some water quality points to return to the equilibrium position is mainly surface tension, which is called capillary waves. In general, the exchange of heat energy between air and seawater is not taken into Consideration. In most cases, the wave spectrum is supposed to be in an ideal state of sufficient growth conditions [5][6][7]. Under this condition, if the wind with constant velocity and constant direction blows across the sea surface for a long time to stabilize the sea, the growth state of sea surface under such an ideal condition is relevant only to the friction wind speed, which is called the steady state spectrum.
In previous research work using numeric simulation of 2-D sea surface based on different ocean wave spectrum models [8], Monte Carlo and fractal method are the most important and most commonly used methods. The 2-D sea surface can reflect the energy distribution characteristics of the wave power spectrum in different directions. The Monte Carlo method is a type of linear filtering method which adopts the methods of fast inverse Fourier transform and power spectral filtering to numerically calculate the 1-D rough function of the ocean wave, thereby realizing the simulation. At the same time, the reconstruction of the sea surface based on the fractal function can well reflect the nonlinear characteristics of the wave, which can not only show the periodic order of the large scale, but also reflect the randomness of the small scale. This method is in line with the physical ocean description method of nonlinear sea surface. In this study, based on the most commonly used methods of Monte Carlo method and fractal method, the linear sea surface and the nonlinear sea surface are simulated and calculated, and the adaptability of the two is preliminarily studied. The JONSWAP spectrum is the most classic spectrum for constructing sea surfaces based on Monte Carlo method [9]. The JONSWAP spectrum is an international standard wave spectrum, which is an unsteady wave spectrum based on statistical analysis of measured data from the Joint North Sea Wave Project. With some shortcomings, however, the spectrum is not applicable to offshore underdeveloped wind and wave spectrum [10]. Fully developed spectrum refers to the state where the wind transfers energy to the wave and the energy dissipated by wave breaking and vortex motion reaches a balance when the fetch and time are large enough at a given wind speed, at which time the wind wave stops growing. The state before reaching this growth state is called an underdeveloped spectrum state. Pierson and Moscowitz conducted spectral analysis experiments on the data of the North Atlantic history, and proposed a general fully developed factorial spectrum, which is called PM spectrum. As a fully developed wave spectrum, the PM spectrum is only suitable for fully developed wind waves, whose high frequency portion is proportional to 5   , where  is the peak shape parameter of the wave spectrum. However, most of observations and Phillips and Kitaĭgorodskiĭ's theoretical results indicate that its high frequency portion should be proportional to 4   . The high-frequency part is also proportional to 5   although the storms at different growth stages can be described by the JONSWAP spectrum, showing inconsistency with most observations and theoretical studies. Moreover, the JONSWAP spectrum limited by the fetch state cannot be used for limited water directional distributions as a deep-water wind wave spectrum. Deep water waves are waving whose wavelength is much smaller than the depth of sea water. After the research of existing wave energy spectra, the theoretical wave power spectrum was deduced by physical oceanographer Wen according to the actual situation of China's offshore wind wave spectra, and the low frequency part was modified by him to enhance the consistency of the spectrum representation and China's offshore sea conditions [11]. The wave spectrum was adopted by the "portal hydrological code" gaining the approval of the ministry of communications of China. The Wen's spectrum is close to the JONSWAP spectrum, but the accuracy of the Wen's spectrum is found to be higher than that of the JONSWAP spectrum after many experiments. Meanwhile, the Wen's spectrum is suitable for different depths and stages of wave growth [12]. According to the inspection of wind and the observation data of waves in most seas of China, the improved Wen's spectrum is more consistent with the measured data compared with other wave spectra. Therefore, in the direction function comparison experiment of this study, based on fifty one groups of data, the adaptability of the Wen's direction spectrum and the Donelan direction function is analyzed. Correspondingly, fractal means the fact that wave components are somewhat similar to the whole and the fractal represents the selfsimilar structure nested inside the shape. Large gravity waves and small capillary waves are considered to constitute the actual sea surface; the characteristics of large-scale order and small-scale disorder can be considered by the fractal model, which makes it feasible to simulate the random surface with the fractal method [13]. Jaggard, Sun and Franceschetti proposed to simulate rough sea surface with a fractal function and constructed the 1-D and 2-D fractal models based on band-limited Weierstrass fractal function [14]. The Weierstrass function approach is applied to investigate the fractal of geometric phase in the wave spectrum. On the basis of the above, the PM power spectrum was proposed by Pierson and Moskowitz according to the spectral analysis of 460 ocean observations in the North Atlantic. First, obtaining the average distribution of random waves, then performing dimensionless and mathematical fitting process on the above distribution values and finally getting the wave number spectrum model of the main wave direction [15], the method is characterized by the fact that a peak can always be seen in the power spectrum regardless of wind speed value, indirectly proving the impact of at least two spectra on the sea surface.
The following is the structure of this study. In the second section, based on the measured wave spectrum distribution data at the Bohai buoy site, there is analysis of the distribution characteristics of measured data, Wen's direction spectrum and Donelan function from the perspectives of directional distribution function and statistical results of the maximum under the conditions of different wind speed and different frequencies of different directional functions. In the third section, a model for fast numeric calculation of 2-D sea surface is proposed based on the above experimental conclusions and in combination with the Donelan function, Wen's power spectrum and Monte Carlo model. In addition, the modified double power law method is applied to the basic theory of Weierstrass self-affine fractal function to make numeric calculation of the 2-D sea surface and to analyze the numeric distribution law of 2-D sea surface in terms of different wind directions, wind speed and fetches. This proves the adaptability of the above-mentioned wind wave calculation model under the conditions of different directional distributions and different wave stages. Finally, by comparing the autocorrelation fitting characteristics of direction distribution under different wind directions and the correlation coefficient of the contour distribution under different wind speeds, there is analysis of the adaptability of the above numeric calculation models under the conditions of different water directional distributions, different sea conditions and different wind and wave growth stages to come to the experimental conclusions. In the fourth section, the contents of above experiments are summarized, the research conclusions are drawn, some problems in the research process are put forward and the future research direction is pointed out.

Experimental Analysis of Distribution Characteristics of Different Directional Functions
As the basic nature of the ocean wave, the direction spectrum is characterized by its time-varying complexity, difficulty in observation and difficulty in indoor experiments. Therefore, there is a lack of marine measured data on direction spectrum and also a lack of experimental verification of related theories. Among the several models of the directional spectrum proposed internationally, most theories use the directional spectrum as the product of the spectrum and the directional function to approximate the directional distribution. That is, is the wave spectrum and ( ) G  is the directional function. Thereby, the empirical approximation process can be realized and the calculation efficiency is greatly improved on the basis of conforming to the physical meaning.
In the field of actual wave direction spectrum research, the accuracy of the JONSWAP spectrum is lower than that of the Wen's spectrum suitable for different depths and wave stages after a large number of experiments [16,17]. According to the inspection of wind wave observation data in most seas of China, improved Wen's spectrum shows higher consistency with measured data than other wave spectrum.
The Wen's direction spectrum is based on analytical methods. We calculate its spectral distribution at different propagation directions and conclude that, ( , ) ( , ) Correspondingly, based on the second-order transformation of the Fourier transform, the Donelan directional spectrum can be obtained from the following formula, In the above formula, / p f f represents the relative frequency and  is the relative azimuth angle. . P is a dimensionless quantity that represents the growth factor and also represents the spectral sharpness, which reflects the growth of the wind wave to some extent. p  represents the main peak frequency, L    represents the lower limit of the angular frequency under the equilibrium domain condition, 0 m represents the energy spectrum at zero order in formula (2), U is the wind speed at 10 m above the surface, g represents the gravitational acceleration and n represents the wave vibration value. Other parameters are a function of relative azimuth and growth factor. The Bohai directional spectrum data measured by the offshore buoy site was selected for analysis. Based on fifty one groups of data, the linear fitting accuracy of the directional spectrum distribution of the Donelan function and the Wen's spectral function in China's offshore waters is analyzed from the perspective of standard deviation of directional distribution and statistical results of the maximum.
In this work, the directional characteristics of different spectra functions are expressed based on the following formula, where the max directional distribution and the standard deviation of the directional distribution function are denoted as respectively. That is, The experimental data are selected from the set average of fifty one groups of data from Bohai 8# buoy site. The measurement method selected during the calculation is called instrument array method. When the data are measured, the external conditions are as follows: the test site water depth is 26.9 m and the array sampling interval is 0.248 s, the standard total length of sampled samples is 19.7 min. Where U represents the average wind speed at 10 m above the surface. In the whole experimental test, the range of the effective wave period In addition, the periodic distribution of the effective wave can be derived by the formula 1 3

1.15
The statistical experimental data are as follows, During the fitting process, standard conversion is needed to facilitate the comparative research.
Combined with the measured data, the directional distribution function , the fitting precision of the different directional distribution function and the measured Bohai buoy data are quantified. The form is as follows.
In the fitting experiment, the error value is calculated by the optimal frequency of the statistical distribution mode. The smaller the error between the two, the higher the fitting accuracy. Combined with different methods, the data distribution adaptability of the measured directional function is summarized, and the following statistical results are obtained. The fitting accuracy of the above directional function and the measured data are as follows [2], wherein, . The smaller the error value, the better the frequency. If there is a minimum number of parallels, the parallel is considered to be optimal. Based on the above principles, the fitting adaptability of the Wen's spectrum and the Donelan function to the measured wave direction spectrum is studied, and the above different distribution patterns are counted, and their errors are calculated. It can be seen from Table 1 that under different sea conditions, when the relative wind speed / p U C is relatively stable, the standard deviation and the highest value of the directional distribution function of the Wen's direction spectrum and the Donelan function under steady-state sea surface decrease with the increase of the relative frequency / p f f . When the relative frequency is 1.0, the above reference value is approximately equal to the lowest value, and then gradually becomes larger, but the increase is weakened, which shows commonality. At the same time, as the relative wind speed gradually increases, the sea state level gradually increases, and the directional distribution function maximum and standard deviation of the Wen's direction spectrum also gradually increase and exhibit anisotropy. On the contrary, the statistical parameters of the Donelan directional function remain stable, are not affected by the sea conditions and have good adaptability. It can be calculated from Tables 1 and 2 that when 0.90 p f f  , the distribution pattern of the Wen's direction spectrum and Donelan function reaches a minimum value and the threshold interval is also minimized when the roughness length is constant. As the wave growth state increases, its distribution range is relatively stable, and the directional distribution function maximum and directional distribution function standard deviation gradually become stable too. From the perspective of wave mechanics, this is actually because the Donelan function characterizes the target spectrum characteristics mainly by referring to the distribution characteristics of the relationship between wave steepness and age. At the same time, the weight function of the wave age state is lower. In contrast, the Wen's spectrum is based on / p U C to represent the rough wave length. The model ignores the influence of whether the sea surface is fully developed to some extent. In fact, the measured data are more accurate and its distribution in different directions is more stable. Moreover, from the actual fitting experiments in Figures 1 and 2, the standard deviation of the directional distribution of the Donelan function is more stable, and the function distribution value is more in line with the theory and its fitting accuracy with the measured data are much better. Therefore, in the experiment process of selecting the direction function to calculate the 2-D sea surface, the Donelan direction function is selected to reflect the anisotropy of the 2-D ocean waves.

Numeric Simulation of Fast 2-D Sea Surface Based on Monte Carlo Method
In At the same time, in the Monte Carlo filtering process, the fundamental amplitude component of the harmonic should be taken as an independent Gaussian variable [18]. In this process, the variance is proportional to power spectrum ( ) j S k . Therefore, the sample set form of the random rough sea surface is expressed as follows [19], In the above formula, the formula for calculating the spectral sharpness factor is the real space for one-to-one correspondence [19]. Fourier's reference equation coefficient satisfy the relationship as follows . After that, using the 2-D IFFT calculation model to implement the formula (6), the numeric experimental results can be obtained. The comprehensive formula for the 2-D sea surface modeling used in above calculation process is as follows, Wherein,  10 U is the wind speed 10 meters from the random sea surface and h is the water depth. The meaning of other parameters is consistent with the parameters in the other formulas above. In the numeric simulation experiment, the sea surface growth state background parameter of the measured 320# data of ice multiparameter imaging X-band radar in 1993 was used as the modeling standard. Among them, sea air humidity, sea surface temperature and downwind direction are 81°,4.3 °C and 320°. The grazing angle is 1.75° and the viewing angle coordinates are (60°, 80°). By adjusting the wind speed to realize the numeric calculation of the wave rough fluctuation of the wind-driven 2-D sea surface, the experimental results are as follows. Based on the above experimental results, it is indicated by Figures 1-3 that there is relatively natural superposition of gravity waves and capillary waves on the sea surface when the wind speed at 10 m from the sea surface is 11 km/h. Therefore, it can be inferred that the sea is at a fully developed state under the fetch condition of 12 km. When the wind speed becomes 16 km/h, the wind fetch remains unchanged at this time, but as the rough surface undulation increases, the rough length of the simulated random rough surface also increases. Relatively speaking, the details of the tension wave and the capillary wave are lost. However, as can be seen from Figure 3, when the wind fetch reaches 15 km, even under the condition of low or medium speed sea state, the undulation of the random surface is large, indicating that the wind fetch is positively related to the rough length of the simulated sea surface.

2-D Sea Surface Computing Model of Weierstrass Self-Affine Fractal Function Combined with Modified Two-Sided Power Law
From the perspective of fractal theory, the random sea surface can be regarded as a combination of macroscopic large-scale structure and subtle small-scale structure. The distribution pattern of the waves is to superimpose small-scale splashes and bubbles on large scale waves of approximately periodicity. This self-similar theoretical model can take into account the large scale gravity waves of random rough waves to maintain orderly features and the small-scale tension waves to exhibit disordered features. Therefore, the fractal function has also become an important theoretical basis for constructing a sea surface model. In the process of studying fractal theory, Mandelbrot proposed that the Weierstrass function has fractal features and its one-dimensional expression is as follows, In the simulation experiment, in order to maintain the scientific rationality of the sea surface numeric calculation results, the background parameter of the sea surface growth state of the 320# measured data of the experimental parameter ice multiparameter imaging X-band radar is also taken as the simulation standard. In order to make the comparison experiment more accurate, the parameters of sea air humidity, sea surface temperature and downwind wind direction are also consistent with the above Monte Carlo modeling method, which are 81°, 4.3 °C and 320°, respectively. At the same time, the average speed of the observed radar is , the time is = 0 s t , the grazing angle is 1.75°, the fractal dimension D is 1.5 and the angle of view is (60°, 80°). The numeric calculation results of the simulated 2-D sea surface at different wind speed based on the above improved classification model are as follows (Figures 4-6).   When the scale factor b is a rational number, the formula (7) is expressed as a standard periodic function. Under this condition, the above sea surface simulation results of Figures 4-6 are characteristic sea surfaces with periodic gravity wave cross. When the scale factor b is an irrational value, the value of the fractal model function becomes a quasi-periodic function. After analyzing the step height value of the wave height, the roughness of the fractal sea surface model can be expressed by the fractal dimension, that is, 2 1 D   . Comparing the simulation models at different wind speed above, the greater the wind speed, the more severe the rough surface changes on the fine structure and the coarser the corresponding rough surface, the higher the height of the sea surface is. Under the condition that the fractal dimension is kept constant, the fluctuation of the random rough sea surface is positively correlated with the wind speeds (Figures 4 and 5) and the fetches (Figures 4 and  6), which is largely similar to the 2-D rough surface simulated by Monte Carlo method, and has certain physical significance. Of course, it should be noted that, limited to the requirements of the length of the study and the setting of the experimental scheme, this study does not study the influence of fractal dimension and wind direction on the height fluctuation of the sea surface.

Comparative Study of Different Numeric Calculation Methods on Sea Surface
In the simulation of the ocean wave spectrum, the inverse spectral Fourier transform is corresponding to the autocorrelation function of the simulated sea surface wave height. Besides, the correlation difference method can also be adopted for more effective study of the distribution difference between the target spectrum and the simulated sea surface spectrum. In this work, we try to measure the correlation between the simulation method and the target spectrum by extracting the correlation difference statistic in the simulated ocean wave spectrum structure corresponding to each 400 × 400 random seed matrix in the simulated sea surface in . This is to analyze the adaptability of different sea spectrum and different numeric simulation methods when simulating the sea surface. First, the wave height corresponding to the highest spectrum is sequentially made to be a difference, and then the average smoothing is performed on the direction distribution of the corresponding wind fetch distribution point, and the obtained set average linear wave is Hilbert transform, that is, Wave height/z residual by determining whether the wave height correlation of two adjacent error terms is 0. This test is based on the assumption that the errors are all generated by the first-order autoregressive process. Figures 7 and 8 show the frequency-domain characteristic curves of the correlation difference distribution in the main wave direction and its orthogonal direction for different sea-field simulation models. Among them, the random phase value is between 0-2π, the wind speed is 11 km/h and 16 km/h, respectively, and the wind zones are 12 km and 15 km, respectively.  Based on the frequency domain characteristics of the correlation difference distribution in the main wave direction and the orthogonal direction of the above different sea-field simulation models, this study analyses the similarity between the simulated sea surface and the target spectrum by analyzing the wave height energy distribution of different numeric calculation methods in the 320° direction. Among them, the default random phase value is between 0-2π, the wind speed is 11 km/h, and the wind fetch is 15 km. At this time, the sea surface is fully developed. The distribution results are as follows (Figures 9 and 10),   The error test for the standard deviation regression model is usually based on the assumption that the i-th error i  and the j-th error j  are irrelevant in the traditional sense, therefore, the correlation in the error term implies that the model established now does not express all the information contained in the data. The reason for the above problem is that the correlation hypothesis test is applied to the modeling of the ocean wave spectrum, and the residuals of adjacent data tend to be similar in time or space. For example, the random phase of the neighborhood of the adjacent wind fetch has a period, which leads to the data tending to have relevant errors, and the error rate also increases regularly with the increase of the simulation step size. The error gain is affected by the common external environment, which is mainly affected by the fetches and wind speeds. It should be specially noted here that the reason why the error gain is not easily affected by the wind direction is that the simulation model itself is established under the ideal sea surface state, and the wind direction remains unchanged at 320°. By comparing the frequency domain characteristic curves of the sea surface correlation difference distribution along the main wave direction and its orthogonal direction in Figures 7 and 8, it is obvious that there is a wider range distribution of the Monte Carlo method than the fractal method in terms of the normalized frequency domain distribution range. This indicates that in the process of simulating the sea surface, when the sea surface state is at an insufficient growth state where the wind speed is lower than 11 km/h and the wind fetch is less than π/6 π/3 π/2 2π/3 5π/6 π 7π/6 4π/3 3π/2 5π/3 11π/6 2π Random phase Normalized energy distribution of wave height π/6 π/3 π/2 2π/3 5π/6 π 7π/6 4π/3 3π/2 5π/3 11π/6 2π Random phase Normalized energy distribution of wave height or equal to 12 km, the sea surface simulated by the Monte method is more prone to rough length fluctuation, while the fractal method is not sensitive. Moreover, the normalized average energy distribution of random phase in Monte Carlo modeling in Figure 9 is more uniform along the length of a period, which is consistent with the anisotropy of the actual wave energy distribution. In addition, the wave height energy distribution curve in Figure 10 is related to the simulated sea surface correlation distribution curve when the wind speed is 11 km/h in Figure 8 and 15 km in the wind area. It also shows from the side that the fractal method still maintains the overall energy part and the rough fluctuation of the wave when simulating the sea surface. The phase is similar, but it does not fully reflect the anisotropy of the gravity wave and the tension wave in the superposition direction of 320°. Therefore, it does not satisfy the sea surface modeling conditions.

Conclusions
Based on the fifty one groups of data on direction distribution measured from buoy sites, the directional distribution characteristic and statistical results of the maximum of Wen's directional spectrum and Donelan function are calculated and analyzed at different relative wind speeds and different relative frequencies in this study. Afterwards, this study makes statistical analysis of the optimal frequency of Donelan function and Wen's directional spectrum under the minimum error value in the process of fitting experiment. Combined with the measured data, this study analyses the numeric fitting characteristics between the measured data and the random rough sea surface with different direction spectra from the perspectives of standard deviation of directional distribution function and statistical results of the maximum under conditions of different wind speed and different frequencies. Thus, it is proved by means of experiments that the directional distribution function standard deviation of the Donelan function is more stable, its function distribution value is more in line with the theory, and the fitting accuracy with the measured data are better. Based on the experimental conclusions, a model for fast numeric calculation of 2-D sea surface is proposed in this work. Then, the study applies the modified double power law to the basic theory of Weierstrass selfaffine fractal function, numerically calculates the 2-D rough surface, and analyses the numeric calculation results of 2-D sea surface under different wind directions, wind speeds and wind fetches. This proves that the above mentioned wind wave calculation model has good adaptability in different water directional distributions and different wave stages. Finally, this study compares the autocorrelation fitting characteristics of the directional distribution of the sea surface simulated by two numeric calculation models under different wind direction angles. At the same time, the correlation coefficient of the contour distribution of the sea surface simulated by two numeric calculation models under different wind speed is compared too. Through the above experiments, this study further analyses the adaptability of the above numeric calculation models under different wind and wave growth stages and points out that the fractal method can keep the overall energy part similar to the rough length of the waves, but it cannot fully reflect the gravity wave and the anisotropy of the tension wave in the superposition direction of 320°.
Moreover, it should be noted that the modeling method based on the Monte Carlo standard spectrum has two distinct advantages as compared with the superposition fractal method. First, the calculation efficiency is relatively high. It can reduce the calculation burden in direct inverse Fourier transform as tree hour can be saved during calculation of the actual sea surface with the wind fetch of 2000 × 2000. Second, the detailed features of sea surface generated by the simulation are more obviously. Because there is a wider distribution range of the Monte Carlo method than the fractal method in terms of the normalized frequency domain. This proves that, in the offshore sea surface simulation, the Monte Carlo method is more suitable for the numeric calculation of the undeveloped sea surface than the fractal function method under the condition of middle and low wind speed and small and medium-sized wind fetch. In addition, the random phase distribution simulated by Monte Carlo method is in the range of 0-2π and the energy distribution in different quadrants is more even, which is consistent with the actual sea surface variation of peaks and troughs alternating with phase in one cycle. In contrast, when the random phase is 0-π, the energy distribution of the fractal function simulation sea surface is all above the zero boundary, which is represented by the sea level fluctuation of the peak and its neighborhood. when the random phase is π-2π, the energy distribution is all below the zero boundary, which is represented by the sea surface fluctuation of the valley and its neighbors. The above state obviously does not match the actual situation. In contrast, the Monte Carlo method modified in this study can keep the overall energy part consistent with the rough length of the waves and can also reflect the anisotropy of gravity waves and capillary waves in the superposition direction.