On the Deterministic Prediction of Water Waves

Abstract

This paper discusses the potential of deterministic wave prediction as one basic module for decision support of offshore operations. Therefore, methods of different complexity—the linear wave solution, the non-linear Schrödinger equation (NLSE) of two different orders and the high-order spectral method (HOSM)—are presented in terms of applicability and limitations of use. For this purpose, irregular sea states with varying parameters are addressed by numerical simulations as well as model tests in the controlled environment of a seakeeping basin. The irregular sea state investigations focuses on JONSWAP spectra with varying wave steepness and enhancement factor. In addition, the influence of the propagation distance as well as the forecast horizon is discussed. For the evaluation of the accuracy of the prediction, the surface similarity parameter is used, allowing an exact, quantitative validation of the results. Based on the results, the pros and cons of the different deterministic wave prediction methods are discussed. In conclusion, this paper shows that the classical NLSE is not applicable for deterministic wave prediction of arbitrary irregular sea states compared to the linear solution. However, the application of the exact linear dispersion operator within the linear dispersive part of the NLSE increased the accuracy of the prediction for small wave steepness significantly. In addition, it is shown that non-linear deterministic wave prediction based on second-order NLSE as well as HOSM leads to a substantial improvement of the prediction quality for moderate and steep irregular wave trains in terms of individual waves and prediction distance, with the HOSM providing a high accuracy over a wider range of applications.

1. Introduction

Ships and offshore structures are exposed to the sea, limiting the scope of application of offshore operations, from efficient and economic offshore operations in moderate sea states to reliability as well as survival in extreme wave conditions. To classify offshore operating conditions, predefined limiting criteria such as absolute or relative motions are typically combined with a limiting characteristic wave height or sea state based on stochastic analysis in the design process, allowing identifying feasible and infeasible sea states, respectively. Typically, these limiting sea states are further reduced in terms of allowable significant wave heights by insurer and surveyor of offshore operations. This approach limits the applicability of offshore operations strictly as most sea states (in particular in the transition area between feasible and infeasible region of a scatter diagram) will feature favourable wave sequences allowing short-term offshore operations which are elapsed unused. Knowing the future in terms of deterministic prediction of the encountering wave field and structure motion a sufficient time span in advance would allow reducing the downtime as the critical part of an offshore operation typically corresponds to a short time interval of the whole operation. In addition, a decision support system (DSS) based on deterministic wave prediction can also detect critical wave groups in pretended feasible sea states increasing the safety of complex operations. The impact of such a DSS on shipping as warning system (e.g., parametric roll, extreme wave events, extreme motions and structure response) should not be left unmentioned even though the focus may lie on offshore operations.

The established research on deterministic wave prediction comprises two indispensable constituents: sea state registration and wave prediction. The deterministic structure motion prediction based on the wave prediction is a further key element for a DSS but will not be discussed in detail within this paper. For the registration of the future encountering wave field, the surface elevation has to be measured at a certain distance from the prediction point (e.g., location of offshore operation). Different measurement methods are available, from point measurements (e.g., wave buoy) in time domain to surface elevation snapshots taken from the ship’s X-band radar in the space domain. As the focus of the paper lies on wave prediction exclusively, the following brief description of the sea state registration methods are only introduced in order to support the approach presented in this paper without any technical details.

Point measurements are disadvantageous compared to a surface elevation snapshot as a single point measurement cannot provide information regarding the directionality of the sea state. This problem may be solved by using a large array of several point measurement devices in order to identify the directionality of the wave components. However, as the availability of point measurement devices in the open ocean are generally limited, this approach would denote to install such an array prior to the offshore operation excluding a plug-and-play application. In addition, this approach is impractical for DSS of ships or offshore structures with forward speed. A further drawback is the principle of measurement as point measurements denote that the wave field is recorded at a fixed location over a specific time. As a consequence, the prediction time (time between the moment at which the prediction is available and the moment at which the predicted waves physically arrive at the prediction location) is reduced by the recording time, independently of the computational time needed for the used wave prediction method [1,2,3,4].

The ship’s X-band radars provide the surface elevation for domains of several square kilometres including information on directionality. Theoretically, the detected wave field can be available after every revolution of the radar antenna (<2 s), even though time is needed for the inversion of the radar clutter (due to the presence of the waves) to the wave field. In contrast to single point measurements, where a sufficient measuring time is needed for an accurate and sufficiently long prediction, this time lag is only limited by the numerical efficiency of the used algorithm. Unfortunately, the principle of measurements of X-band radars also yield a main drawback as the radar cannot detect wave troughs behind (steep) wave crests. Thus, the accuracy is dependent on the wave inversion algorithm. However, the advantages of using a ship’s X-band radar (plug-and-play, very fast wave registration over a large space domain, directionality of the wave field, and large prediction horizon) seems to outweigh the disadvantages as two commercial prediction systems are available from the companies Applied Physical Sciences Corporation (FutureWaves™ [5]) and Next Ocean™ [6]. Even this paper focusses on wave prediction exclusively, the presented research assumes surface elevation snapshots (taken continuously by a ship board radar at great distance ahead the operational area) as input for the deterministic wave prediction.

One principal point of the development process of deterministic wave prediction is the implementation of very fast algorithms to obtain the predicted wave field in a sufficient time span before the physical wave field arrives at the prediction location. Only a few methods are capable due to the contrary specifications of very fast calculation time and high accuracy at once. The fastest method but also simplest one is the linear theory, which has already been applied for wave prediction applications with promising findings [2,7,8,9,10,11,12,13,14,15]. A Shipboard Routing Assistance system (SRA) based on the continuous ship’s X-band radar measurements of the surrounding seaway were presented by Payer and Rathje [7]. The linear evolution of continuously measured surface elevation snapshots using the ship’s X-band radar are also the basis for a decision support system for Computer Aided Ship Handling (CASH) developed by Clauss et al. [8]. It was shown that this tool can predict the encountering wave field as well as the structure response fairly accurate for moderate long crested sea states. Later, Clauss et al. [12] and Kosleck [2] enhanced this tool to predict the short crested encountering wave field and the corresponding 6DOF motion. Tyson and Thornhill [16,17] applied linear wave prediction using WaMoS II™ on full scale tests evaluating the predictive skills of the linear methods. They showed that for short forecast duration and small to moderate wave steepness, the accuracy of the linear approach is sufficient. In addition, Naaijen and Huijsmans [9,11] as well as Naaijen et al. [10,15] applied linear wave evolution equations for real time wave prediction and ship motion estimation in long as well as short crested waves concluding “that a 60 s accurate forecast of wave elevation is very well feasible for all considered wave conditions and motion predictions are even more accurate”. Due to the simpleness and robustness, the linear approach is an integral part of commercially available prediction system (e.g., FutureWaves™ [13,14], Next Ocean™ [6]).

However, the linear approach implies uncertainties due to its strong simplifications of the water wave problem. Non-linear effects become dominant with increasing wave steepness reducing the accuracy of the linear approach significantly. Non-linear methods enable advanced simulations with high accuracy but at the expense of computation time. One of the few non-linear and numerical efficient methods suitable for wave prediction are the so-called envelope equations of the NLSE framework [18,19,20,21]. Both the NLSE and the modified NLSE (MNLSE) were extensively experimentally and numerically investigated in terms of non-linear wave evolution (e.g., [22,23,24,25,26] among others). Even if the NLSE captures relevant non-linear phenomena and provides good accuracy for sufficiently narrow spectra, it was generally shown that the NLSE is less suitable for the prediction of irregular sea states due to the spectral bandwidth constraint. The MNLSE on the other side provides significant better agreement due to the broadening of the bandwidth constraints (e.g., [22,23,24,25,26]).

Ruban [27] showed that the NLSE is suitable for the prediction of extreme waves applying the Gaussian variational ansatz to the NLSE in order to obtain a semi-quantitative prediction of non-linear spatio-temporal focussing. Farazmand and Sapsis [28] presented a reduced-order prediction of rogue waves using MNLSE. Their procedure is based on the identification of elementary wave groups (EWG) assuming that critical EWGs do not interact with other EWGs within the short-term prediction. The measured irregular sea state is decomposed into EWGs and the evolution and amplitude growth of critical EWGs is determined by precomputed EWGs based on MNLSE. Thus, the numerical effort for the prediction of extreme events is reduced to the proper decomposition of the sea state and identification of the underlying, precomputed EWG. Cousins et al. [29] showed a similar approach: a data-driven prediction scheme for the prediction on extreme events. This approach can also divided into two separate components. The first component provides the data base for the prediction of extreme events, is based on the MNLSE and can be applied prior any real world application. The basis are localized wave groups which are numerically evolved for varying amplitudes and periods resulting in specific group amplification factors. The second component is the real-time extreme wave predictor using field measurements. Within the measurements, coherent wave groups are to be identified and the future elevation is estimated from the data base. Again, the numerical effort is reduced to decomposition of the irregular wave sequence, identification of coherent wave groups, and resulting amplification.

The NLSE as well as MNLSE has also been used for deterministic wave prediction (e.g., [26,30]). Simanesew et al. [30] applied both for the deterministic prediction of long- as well as short-crested sea states in time domain. They showed that both the NLS and MLSN provide sufficient accuracies in long-crested sea states for short forecast distances and horizons. For increasing directionality the quality of the forecast decreases significantly. Of particular note in this investigation is that the linear dispersion behaviour is captured exactly by using the dispersion operator introduced by Trulsen et al. ([21]).

An alternative non-linear method for wave prediction is the numerically efficient high-order spectral method (HOSM) [31,32,33]. Wu [1] as well as Blondel et al. [34] applied the HOSM for the deterministic reconstruction and prediction of non-linear wave fields. The basis of their work are wave registrations based on one or several wave probes at specific locations. Thus, sophisticated optimization procedures for the reconstruction of the wave field in space were used at the beginning as snapshots in the space domain are required as initial condition for the HOSM. This time-consuming reconstruction process hindered an effective application of the HOSM for wave prediction. However, taking advantage of the supposed drawback by using surface elevation snapshots continuously taken by a ship’s X-band radar at great distance as input, enables the use of the efficient and accurate HOSM without time-consuming reconstruction of the wave field in space domain. Nevertheless, they generally showed that long-time and large-space simulations of non-linear sea state evolutions can be performed accurately and efficiently with the HOSM. Clauss et al. [35] and Klein et al. [36] presented that the HOSM predicts non-linear wave group propagation very accurately. The applicability of the HOSM for non-linear real-time prediction was shown by Köllisch et al. [4]. Desmar et al. [37] applied the HOSM for the generation of reference wave snapshots in order to evaluate reconstruction methods of different complexity and the consequences on the prediction accuracy. This investigation showed that non-linear methods for wave reconstruction (based on spatio-temporal optical measurements) are improving the accuracy of the reconstructed initial wave field and prediction.

This paper presents a comparative study on the accuracy of intended wave prediction methods of different complexity. The objective is the evaluation of the applicability of the utilised methods for accurate deterministic wave prediction using irregular sea states. The focus lies on the investigation of the influence of the wave steepness, wave propagation distance and shape of the underlying spectrum. To be consistent with a possible real world application (input snapshots from ship’s X-band radar), the numerical and experimental investigations are performed in space domain. Herein, the unavoidable radar shadow around the ship’s position is also modelled and the influence is evaluated. For the experiment’s validation, a semi-experimental validation procedure is introduced as measuring the surface elevation in space domain is almost impractical in the controlled environment of a seakeeping basin. The accuracy of the predictions is evaluated with the surface similarity parameter allowing an exact, quantitative evaluation.

2. Wave Theory

We now present a brief description of the water wave problem without attempting to be comprehensive, herein only those methods are presented which are relevant for this paper. At the beginning the boundary value problem of potential flow theory is briefly introduced, followed by the used methods presented in the order of complexity, from linear wave theory via weakly non-linear envelope equations to (fully) non-linear simulations. At the end, the fully non-linear numerical wave tank waveTUB is additionally presented as waveTUB depict the key element for the semi-experimental validation procedure. For this study, the water wave problem is simplified to long-crested waves.

Assuming that the Newtonian fluid is incompressible, inviscid and irrotational, the evolution of long-crested waves (travelling in x-direction) can be mathematically described by the following governing equations:

$$▵\mathsf{\Phi }=0,$$

(1)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_z-{\mathsf{\Phi }}_x{\zeta }_x-{\zeta }_t=0& \mathrm{on}z=\zeta (x,t),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill gz+\frac{1}{2}{\left(▿\mathsf{\Phi }\right)}^2+{\mathsf{\Phi }}_t=0& \mathrm{on}z=\zeta (x,t),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_z=0& \mathrm{on}z=-d,\hfill \end{array}$$

(4)

with $▵\equiv ({\partial }^2/\partial x^2,{\partial }^2/\partial z^2)$, $▿\equiv (\partial /\partial x,\partial /\partial z)$ and the subscripts represent the corresponding derivation. The bottom and the unknown free surface representing the boundaries of the fluid domain. Equation (4) presents the boundary condition at the bottom which is considered to be level, rigid and impermeable. For the dynamic boundary condition at the free surface (Equation (3)), the dynamic pressure is defined to be constant equalling the atmospheric pressure. The kinematic boundary condition at the free surface (Equation (2)) defines that particles, being element of the free surface at rest, must not leave the free surface in the presence of waves. The kinematic as well as dynamic boundary condition have to be fulfilled at the unknown free surface $z=\zeta (x,t)$ complicating the solution of the boundary value problem significantly.

2.1. Linear Wave Theory

Linear wave theory can be deduced from the Cauchy problem (Equations (1)–(4)) by applying perturbation theory for the unknown potential and surface elevation. In addition, the boundary value problem at the unknown surface can be approximate by Taylor series expansion. For Stokes wave theory, the perturbation parameter is related to the wave steepness $ϵ={\zeta }_{a}k$. Truncation of the series expansion at order $\mathcal{O}\left({ϵ}^1\right)$ (assuming that the wave height is significantly smaller compared to the wave lengths) results in a simplified boundary value problem for which the analytical linear solution can be determined. For irregular surface elevations of long-crested waves in time and space, the solution reads

$$\zeta (x,t)=\sum _n{\zeta }_{an}cos(k_{n}x-{\omega }_{n}t+{\varphi }_n),$$

(5)

with ${\zeta }_{an}$, $k_n$, ${\omega }_n$ and ${\varphi }_n$ as amplitude, wave number, angular frequency and phase of the $n^{th}$ component wave. The desired sea state is regarded as superposition of independent harmonic component waves, each with a particular amplitude, frequency and phase. Applying Equation (5), a measured irregular surface elevation can be transformed at any position in time and space enabling a very fast and simple approach for deterministic wave prediction. The analytical basis for irregular sea states is the Fourier transform. The angular frequency and wave number are coupled via the dispersion relation

$${\omega }_n=\sqrt{k_{n}gtanh\left(k_{n}d\right)},$$

(6)

with water depth d and gravitational acceleration g. The linear description of the natural seaway enables a simple handling of a complex process with widely acceptable results for engineering applications.

2.2. Non-Linear Schrödinger-Type Equations

The classical NLSE can be derived similarly to the linear theory by applying perturbation and Taylor series expansions on the Cauchy problem [18,19]. Assuming small amplitude waves and a narrow bandwidth spectrum, the perturbation parameters are the wave steepness $ϵ=k_c{\zeta }_a$ and the relative bandwidth $\mu =\Delta k/k_c<<1$ (assuming ${\mu }^2\propto ϵ$). The Taylor series expansions about still water level is introduced in order to simplify the surface boundary conditions. Furthermore, the boundary value problem is separated into the different orders $\mathcal{O}\left({ϵ}^{\left(n\right)}\right)$ and each order is solved sequentially from the lowest to the highest order. Truncation of the perturbation expansion at order $\mathcal{O}\left({ϵ}^3\right)$ leads to the NLSE. The NLSE captures some important aspects of non-linearity of water waves such as modulation instability and provides exact solutions enabling fundamental research on weakly non-linear wave-wave interaction. However, the NLSE is limited in non-linearity magnitude and spectral bandwidth hindering a diversified application as deterministic wave prediction method. These limitations can be reduced by including higher order terms of the asymptotic series expansion. The next higher order $\mathcal{O}\left({ϵ}^4\right)$ leads to the well known MNLSE [20,21]. Brinch-Nielsen and Johnson [38] derived the MNLS for arbitrary water depth. However, in most of the cases, the MNLSE is applied assuming infinite water depth and is typically referred as Dysthe equation. Recently, Slunyaev and Pelinovsky [39] extended the Dysthe equation to the next order ($\mathcal{O}\left({ϵ}^5\right)$) referring to higher-order Dysthe equation.

Generally, and not only with regard to this study, the water depth influence cannot be neglected for deterministic wave prediction as offshore operations are conducted in arbitrary water depth. For the NLSE, the derivation of the necessary coefficients taking the water depth influence into account can be found in Mei [40]. For the next higher order, Sedletsky [41] derived a second-order NLSE ($\mathcal{O}\left({ϵ}^4\right)$) for arbitrary water depth, which reduces to the classical Dysthe equation in the limit of infinitely deep water. Based on this outcome, Slunyaev [42] derived the third-order NLSE ($\mathcal{O}\left({ϵ}^5\right)$) for arbitrary water depth investigating the impact of the different non-linearity orders on the modulation instability. For this study on deterministic wave prediction, the NLSE up to the second-order is used in the following ways.

According to Slunyaev [42], truncation of the perturbation expansion at order $\mathcal{O}\left({ϵ}^4\right)$ leads to the temporal evolution equation

$$i\left(\frac{\partial A}{\partial t}+V\frac{\partial A}{\partial x}\right)+{\beta }_1\frac{{\partial }^{2}A}{\partial x^2}+{\alpha }_1{\left|A\right|}^{2}A+i{\beta }_2\frac{{\partial }^{3}A}{\partial x^3}+i{\alpha }_{21}{\left|A\right|}^2\frac{\partial A}{\partial x}+i{\alpha }_{22}{A}^2\frac{\partial A^{*}}{\partial x}=0,$$

(7)

with A being the envelope of the wave sequence, V the group velocity of the carrier wave, the superscript * denotes conjugate complex. The indices of the coefficients $\alpha$ and $\beta$ mark the specific order of the underlying harmonic (first index) as well as order of expansion (second harmonic), which is adopted from Slunyaev [42] in order to be consistent with the presented equations for the coefficients (summarized in Appendix A). The group velocity depends on the carrier wave number $k_c$, the corresponding angular frequency ${\omega }_c=\sqrt{k_{c}gtanh\left(k_{c}d\right))}$ and the water depth d,

$$V=\frac{{\omega }_c}{2k_c}\left(1+\frac{2k_{c}d}{sinh\left(2k_{c}d\right)}\right).$$

(8)

The first four terms in Equation (7) display the classical first order NLSE and the underlying wave field can be reconstructed as follows:

$$\zeta =Re\left(Aexp(i{\omega }_{c}t-i{k}_{c}x)\right).$$

(9)

Truncation at order $\mathcal{O}\left({ϵ}^4\right)$ results in Equation (7) representing the solution for the second order NLSE (${\mathrm{NLSE}}_2$) equivalent to the MNLSE (and Dysthe equation for $k_{c}d\to \infty$). The reconstruction formula reads

$$\zeta =r_{01}{\left|A\right|}^2+Re\left(Aexp(i{\omega }_{c}t-i{k}_{c}x))+Re(r_{21}{A}^{2}exp(2i{\omega }_{c}t-2i{k}_{c}x)\right),$$

(10)

with contributions of the zeroth, first and second harmonics.

The dispersive behaviour of the waves is represented by the second term and all terms related to ${\beta }_{nm}$; the non-linear wave-wave interaction corresponds to the ${\alpha }_{nm}$ terms. The dispersive terms in Equation (7) are the consequence of the Taylor series expansion of the linear dispersion relation so that the dispersion of the waves is represented in the vicinity of the underlying carrier wave number $k_c$. For increasing order, the obtained dispersion behaviour approximates to the full linear dispersion. However, Trulsen et al. [21] introduced the full linear dispersion operator

$$\widehat{L}=\sqrt{(k-k_c)gtanh\left((k-k_c)d\right))-{\omega }_c},$$

(11)

with k representing an array of wave numbers and the higher order dispersive terms can be calculated by

$${\widehat{F}}^{-1}\left\{\widehat{L}\widehat{F}\left\{A\right\}\right\}=iV\frac{\partial A}{\partial x}+{\beta }_1\frac{{\partial }^{2}A}{\partial x^2}+i{\beta }_2\frac{{\partial }^{3}A}{\partial x^3}+{\beta }_3\frac{{\partial }^{4}A}{\partial x^4}+\dots ,$$

(12)

with $\widehat{F}$ marks the Fourier transform and ${\widehat{F}}^{-1}$ the inverse Fourier transform, i.e., the operation is done in Fourier space by means of pseudo-spectral approach. Equation (12) takes terms beyond the highest order ($\mathcal{O}\left({ϵ}^4\right)$) considered in this paper into account and is formally not consistent with the series expansions truncated at specific orders. Nevertheless, introducing Equation (12) into Equation (7) provides the advantage that the linear dispersive behaviour is fully taken into account resulting in increasing accuracy compare to the classical Taylor series expansion of the linear dispersion relation at each selected order. Another advantage is more practical, as in particular for the determination of higher order derivatives via pseudo-spectral methods, the difficulties with numerical noise and aliasing effects are reduced.

For this study, the two non-linear envelope equations resulting from Equation (7) are implemented in order to investigate the influence of the different orders on the wave prediction accuracy. For each order, the full dispersion operator (Equation (12)) is applied. In addition, the Taylor series expansion of the linear dispersion is also applied for the NLSE (first four terms in Equation (7)) to investigate and evaluate the benefit of applying the full dispersion operator. The classical NLSE simulations including Taylor series expansion of the linear dispersion are performed by implementing the pseudo-spectral split-step method. Hereby, the linear and non-linear part of the NLSE are determined separately. The linear part is calculated in frequency domain applying Fourier transform and the non-linear part in time domain. For all other simulations including the full dispersion operator, the linear and non-linear parts are determined together with a pseudo-spectral approach at each time step and advanced in time with the midpoint finite-difference approximation. The carrier ware number $k_c$ necessary to run the simulations is determined for each snapshot by applying

$$k_c=\frac{\int k·E{\left(k\right)}^{P}dk}{\int E{\left(k\right)}^{P}dk},$$

(13)

with E the variance spectrum of the snapshot and $P=5$ as weighting factor to weight the result toward the visual peak [43].

2.3. High-Order Spectral Method

The HOS method was introduced independently by West et al. [31] and Dommermuth and Yue [32] (cf. Tanaka [33]). This procedure enables the non-linear simulation of short-crested sea states and takes all non-linear interactions, resonant and non-resonant, into account. In addition, wave-current as well as wave-bottom interactions can be considered. For our investigation, the numerical procedure presented in West et al. [31] is implemented using a pseudo-spectral method. Specifically, all derivatives related to the potential and surface elevation are determined in Fourier space assuming periodic boundary conditions while the non-linear products are calculated in physical space. “The (near) linear computational effort” as well as “exponential convergence (...) are notable characteristics of the computational efficacy of HOS methods” [44].

At the beginning, the Cauchy problem is converted to equations at the free surface $\mathsf{\Psi }(x,t)\equiv \mathsf{\Phi }(x,\zeta (x,t),t)$. Using chain rules, the free surface boundary conditions (Equations (2) and (3)) can be rewritten as:

$$\begin{array}{cc}\hfill {\zeta }_t=-{\mathsf{\Psi }}_x{\zeta }_x+W(1+{({\zeta }_x)}^2)& \mathrm{on}z=\zeta (x,t),\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_t=-g\zeta -\frac{1}{2}{({\mathsf{\Psi }}_x)}^2+\frac{1}{2}{W}^2(1+{({\zeta }_x)}^2)& \mathrm{on}z=\zeta (x,t),\hfill \end{array}$$

(15)

with $W={\mathsf{\Phi }}_z{|}_{z=\zeta }$ as vertical velocity at the free surface.

Using Equations (14) and (15) as free surface boundary conditions, the boundary value problem is now exclusively related to the vertical velocity $W(x,\zeta (x,t),t)$ (in space domain) and the solution for $W(x,\zeta (x,t),t)$ in terms of $\zeta (x,t)$ and $\mathsf{\Psi }(x,t)$ can be determined by series expansion. The procedure proposed by West et al. [31] starts from the formal expression that the velocity potential $\mathsf{\Psi }(x,t)$ and vertical velocity $W(x,\zeta (x,t),t)$ can be represented as Taylor series expansion a $z=0$. Assuming that $\mathsf{\Psi }(x,t)$ and $\zeta (x,t)$ are quantities of $\mathcal{O}\left({\zeta }^n\right)$, $\varphi (x,t)$ and $W(x,t)$ are expanded by perturbation series, with $\zeta$ as ordering parameter and $M=m+1$ is the order of approximation of non-linearity. Separating the terms of each order $\mathcal{O}\left({\zeta }^n\right)$ yields,

$${\varphi }^{\left(0\right)}(x,t)=\mathsf{\Psi }(x,t),$$

(16)

for the first order. The next higher order solutions are obtained successively from the next lower order solution,

$${\varphi }^{\left(m\right)}(x,t)=-\sum _{n=1}^m\frac{{\zeta }^n}{n!}\frac{{\partial }^n}{\partial z^n}{\varphi }^{(m-n)}(x,t).$$

(17)

The vertical velocity W at the surface $\zeta$ is obtain by

$$W^{\left(m\right)}(x,t)=-\sum _{n=0}^m\frac{{\zeta }^n}{n!}\frac{{\partial }^{n+1}}{\partial z^{n+1}}{\varphi }^{(m-n)}(x,t).$$

(18)

Assuming periodic boundary conditions, the Fast Fourier Transform (FFT) can be used for determining ${\varphi }^{\left(m\right)}(x,t)$ and its derivatives as well as ${\zeta }_x$. For arbitrary water depth, the wave-bottom interaction on ${\varphi }^{\left(m\right)}(x,t)$ can be derived in Fourier space by

$${\varphi }^{\left(m\right)}(x,t)=\frac{1}{2\pi }\underset{-\infty }{\overset{\infty }{\int }}{\widehat{\varphi }}^{\left(m\right)}(k,t)cosh\left(\left|k\right|(z+h)\right)e^{ikx}dk,$$

(19)

with ${\widehat{\varphi }}^{\left(m\right)}(k,t)$ as Fourier coefficient.

For this study, a low pass filter according to West et al. [31] is implemented to avoid Fourier space aliasing for the higher order terms and the series expansion was expanded up to the fourth order ($M=4$). To suppress high frequency contamination occurring for the highest waves, which can cause numerical instabilities [35], an exponential damping term was additionally introduced in the Fourier space. However, the procedure is neither capable of simulating steep waves close to breaking nor handling wave breaking effects. In this study, the fourth order Runge-Kutta-Gill method is applied to advance the evolution equations in time.

2.4. WaveTUB

This potential theory solver was developed at the Technische Universität Berlin (TUB) for the simulation of non-linear wave propagation [45]. The two-dimensional, non-linear free surface flow problem (Equations (1)–(4)) is solved in time domain. The velocity potential is calculated in the entire fluid domain with the Finite Element Method and at each time step a new boundary-fitted mesh is created. Based on the calculated velocity potential, the velocities at the free surface are determined by second-order differences. The fourth-order Runge-Kutta formula is applied to develop the solution in time domain. At one side of the numerical wave tank, a numerical beach is implemented by adding artificial damping terms to the kinematic and dynamic free surface boundary condition enabling long term simulations. On the other side, a moving wall is implemented for the generation of waves enabling the simulation of piston-type, flap-type and double flap-type wave boards. A detailed description of the theoretical background can be found in Steinhagen [45]. The established program waveTUB showed a high accuracy over recent years for a multitude of different tasks (cf. [45,46,47]).

3. Experimental Program and Results

The following study comprises investigation of irregular sea states based on JONSWAP spectra evaluating the influence of wave steepness, the shape of the underlying spectrum (in terms of enhancement factor $\gamma$) and the wave propagation over large distances. The above introduced intended wave prediction methods are used and the accuracy is evaluated quantitatively. A realistic real world application is simulated by taking eight consecutive snapshots per sea state (similar to a ship’s radar) into account. In addition, the unavoidable radar shadow is also modelled and its influence on the prediction accuracy is discussed. Following, all data are presented in full scale (model scale 1:75).

Table 1. Overview of the investigated irregular sea states.

No.	${\mathit{T}}_{\mathit{p}}$	${\mathit{k}}_{\mathit{p}}\mathit{d}$	${\mathit{H}}_{\mathit{s}}$	$\mathit{ϵ}$	$\mathit{\gamma }$
1	$8.58$ s	$4.1$	$1.83$ m	$0.05$	1
2					3
3					6
4			$3.66$ m	$0.1$	1
5					3
6					6
7			$5.48$ m	$0.15$	1
8					3
9					6

presents the investigated sea states based on JONSWAP spectra with varying significant wave height $H_s$ and enhancement factor $\gamma$. Therefore, three different significant wave heights are investigated based on preselected wave steepness $ϵ=k_p·H_s/2$ marking small, moderate and steep irregular sea states. For each wave steepness, the enhancement factor $\gamma$ is varied. A special feature of this investigation is the implementation of an identical phase distribution (at the wave board) for all investigated sea states. This allows an exclusive evaluation of the influence of wave steepness and bandwidth of the underlying spectrum. Thus, differences in the evolution of the sea states (and thus to the accuracy of the intended prediction methods) are exclusively related to the underlying wave spectrum allowing conclusions on the effect of non-linear wave propagation due to wave steepness as well as band width of the spectrum.

The experiments are performed in the seakeeping basin of the Ocean Engineering Division at TUB. The basin is 110 m long and 8 m wide with a measuring range of 90 m. The water depth is 1 m. The waves are generated by an fully computer controlled, electrically driven wave generator, which can be used in piston as well as flap-type mode. The implemented wave generation software enables the generation of regular waves, transient wave packages, irregular sea states as well as tailored (critical) wave sequences. On the opposite side, a wave damping slope is installed to suppress disturbing wave reflections.

The accuracy of the different wave prediction methods applied on the measured irregular sea states is investigated with an semi-experimental approach. The aforementioned real world application—deterministic wave prediction based on ship’s X-band radar—serves as prototype. Thus, surface elevation snapshots in space domain are used as input. In addition, fixed prediction locations are defined at which the prediction in time domain is presented. The influence of the radar shadow on the prediction accuracy is additionally investigated. To avoid time-consuming and expensive series of measurements, which would be necessary to detect the surface elevation in space domain in the seakeeping basin, a semi-experimental approach is introduced for the preparation of the input snapshots. Measuring the surface elevation in space domain would denote that plenty of successive, pointwise measurements have to be conducted (with suitable calm down time in between) to obtain the snapshot.

The semi-experimental approach enables significantly fewer measurements by more validation scenarios at the same time. The fully non-linear numerical wave tank waveTUB is applied for the determination of the required surface elevation snapshots for the deterministic wave prediction. The seakeeping basin, used for the experimental validation, can be modelled exactly with this numerical wave tank (including the wave maker) resulting in identical boundary conditions at the wave board, i.e., providing an identical starting point of the wave evolution for both tanks. Thus, the irregular sea states are investigated in the seakeeping basin and reproduced with the numerical wave tank waveTUB providing the input wave snapshots for the deterministic wave prediction. The input snapshots based on waveTUB as well as the prediction results are experimentally validated.

Figure 1 presents the general scheme of this procedure showing the side view on the wave tank (all data full scale). On the left hand side the wave maker and on the right hand side a wave damping beach is installed/modelled. A surface elevation snapshot at a specific time is also indicated which is the basis for the intended wave prediction methods. For the pointwise validation of the waveTUB input snapshots as well as wave predictions by measurements in time domain at fixed positions, four wave gauges are installed during the experiments. The waveTUB reproductions are compared with the experiments near the wave maker at wave gauge 1 to draw conclusions on the quality of the reproduction. The predicted wave sequences for three different prediction distances are validated at the other three wave gauges (2-4).

The waveTUB simulations are used to extract eight consecutive snapshots from each sea state simulating consecutive snapshots taken from the ship’s radar of a real world application. Each consecutive snapshot is used as input for the wave prediction for the three location (cf. Figure 1 wave gauges (2-4)). Figure 2 presents exemplary the applied procedure in detail for sea state 1. The surface elevation snapshots (grey curves) taken from the waveTUB simulations are presented in the left diagrams—the consecutive time steps are displayed from top to bottom. In dependence on the real world application, the real world radar shadow close to an offshore structure is simulated by modifying the waveTUB results. Assuming that wave gauge 2 is the location of the radar and a realistic radar shadow of 500 m, the radar shadow is simulated by changing the surface elevation to zeros in this area (as well as behind wave gauge 2). The modified surface elevation snapshots are shown as red curves in the left diagrams of Figure 2. In addition, the positions of the four wave gauges used in the experiments are indicated by the black vertical lines.

The measured surface elevations (blue curves) in time domain at wave gauge 2 are exemplary compared with the HOSM simulations (red curves) on the right hand side of Figure 2—the prediction basis are the modified waveTUB input snapshots (red curves on the left hand side). Only the time span where an accurate prediction can be accomplished theoretically is presented. This prediction region marks the spatio-temporal domain which can be predicted based on known ranges in space and time for certain sea states. The group velocity of the fastest and slowest wave groups within the sea state define the prediction range—the slowest wave group must have reached the target location (position of offshore structure or in our study the positions of wave gauges 2-4) and the fastest wave group must not have passed the target location (cf. [1,2,3,48]). For this study, the prediction region is determined for each consecutive input snapshot separately but the variation for the different consecutive snapshot as well as sea states is marginal. Based on this, the average prediction time at wave gauge 2 is approximately $t_{pt}\approx 300$ s (cf. red curves in Figure 2 on the right hand side) with an average minimum forecast horizon time $t_{minfh}\approx 90$ s (time the slowest wave group within the sea state needs to reach the target location, i.e., starting point of the prediction range) and an average maximum forecast horizon $t_{maxfh}\approx 390$ s neglecting any computational time (wave gauge 3: $t_{pt}\approx 257$ s, $t_{minfh}\approx 223$ s and $t_{maxfh}\approx 480$ s; wave gauge 4: $t_{pt}\approx 213$ s, $t_{minfh}\approx 355$ s and $t_{maxfh}\approx 568$ s). Please note that these values depend strongly on the dominant wave length and the underlying wave spectrum, i.e., with increasing wave length the forecast horizon will decrease.

The measured and calculated wave sequences at specific locations in time domain are evaluated quantitatively applying the surface similarity parameter (SSP) [49]. The SSP represents “a quantitative method to compare temporal or spatial series in one dimension or temporal or spatial surfaces in two dimensions” [49]. The SSP is a normalized error between two signals or surfaces written in terms of Sobolev norms

$$SSP=\frac{(\int |{\widehat{F}}_{f1}\left(k\right)-{\widehat{F}}_{f2}\left(k\right){{|}^{2}dk)}^{1/2}}{\int |{\widehat{F}}_{f1}{\left(k\right)|}^2{dk)}^{1/2}+\int |{\widehat{F}}_{f2}\left(k\right){{|}^{2}dk)}^{1/2}},$$

(20)

with $\widehat{F}\left(k\right)$ as Fourier transform of the two signals. The magnitude of the SSP varies between 0 (perfect agreement) and 1 (perfect disagreement). One important advantage compared to other available coefficients is the inclusion of both the amplitude and the phase difference of two series or surfaces [49].

Figure 3 and Figure 4 present the accuracy of the predicted wave sequences compared to the measurements for all investigated methods. For the sake of clarity, the results are divided into Figure 3 comparing all NLS-related findings and Figure 4 comparing the accuracy of linear transformation, second-order NLSE and HOSM.

Figure 3 presents the SSP of the waveTUB input snapshot (black plus), the classical first-order NLSE (orange curves), the first-order NLSE including full dispersion operator (NLSE_FD, magenta curves) as well as the second-order NLSE (NLSE₂, green curves) simulations for the investigated sea states. The darker illustrated curves represent the prediction accuracy including the radar shadow and the lighter illustrated curves display the accuracy without radar shadow. The dots on the curves illustrate the positions of the three wave gauges. The SSP for the predicted sea states represents the mean value of the eight consecutive forecasts. The corresponding variance for the prediction with radar shadow is presented by the transparent areas around the respective curves (with the same colour). The SSP for the waveTUB input snapshot is calculated for the whole time trace at wave gauge 1. Each row presents sea states with constant wave steepness but increasing enhancement factor from left to right and each column sea states with constant enhancement factor but increasing wave steepness from top to bottom. The black vertical line marks the border between “radar” input snapshot (left) and the prediction zone (right).

Noteworthy at a first glance is the influence of the inclusion of the full dispersion operator within the first-order NLSE: the accuracy of the first-order NLSE is significantly improved by taking the full dispersion into account. Based on the underlying simplified assumptions of the NLSE, two main reasons can be identified. On the one hand, the assumed narrow bandwidth of the wave spectrum is not fulfilled for the investigated sea states. Thus, the assumption that the waves within the spectrum evolve with the same group velocity based on a dominant wave length is too inaccurate for deterministic wave prediction of arbitrary irregular sea states. This is clearly supported by the data shown in Figure 3 top as the accuracy of the NLSE increases with increasing enhancement factor (narrower spectrum bandwidth). On the other hand, the range of validity of the NLSE are small amplitude waves. As a consequence, the gain in accuracy with increasing enhancement factor is less distinctive for the steepest irregular sea states. NLSE_FD and NLSE₂ have the same accuracy for the smallest steepness. This shows, that the NLSE can be applied for broader wave spectra by taking the full dispersion into account. With increasing wave steepness, the advantage of the NLSE₂ becomes clear. The accuracy of the NLSE_FD decreases with increasing steepness compared to the waveTUB input snapshot accuracy, whereas the accuracy of the NLSE₂ is indeed also effected but not so distinct showing a better overall performance also for the steepest cases.

Figure 4 presents the SSP of the waveTUB input snapshot (black plus), the HOSM (red curves), NLSE₂ (green curves) as well as linear (blue curves) simulations for the investigated sea states. The diagram is arranged as the previous one, the darker illustrated curves represent the prediction accuracy including the radar shadow and the lighter illustrated curves display the accuracy without radar shadow. The dots on the curves illustrate the positions of the three wave gauges. The SSP for the predicted sea states represents the mean value of the eight consecutive forecasts. The corresponding variance for the prediction with radar shadow is presented by the transparent areas around the respective curves (with the same colour). The SSP for the waveTUB input snapshot is calculated for the whole time trace at wave gauge 1. Each row presents sea states with constant wave steepness but increasing enhancement factor from left to right and each column sea states with constant enhancement factor but increasing wave steepness from top to bottom. The black vertical line marks the border between “radar” input snapshot (left) and the prediction zone (right).

The first row (same small wave steepness and increasing enhancement factor) shows that for small wave steepness the HOSM, the NLSE₂ and linear results are in the same accuracy range. This indicates that applying linear transformation for wave prediction of sea states with small wave steepness is sufficient and there is no need for more sophisticated methods. In addition, the three top diagrams reveal that the accuracy of the prediction remains almost constant over a very large distance, i.e., 2025 m or 18 dominant wave length from the beginning of the input snapshot and 4875 m or 43 dominant wave length from the end of the input snapshot to wave gauge 4. The main reason is the fact that for sea states with small wave steepness non-linear effects are very small and these non-linear effects evolve only in a large time and space scale. However, comparing the diagrams of each column shows clearly that with increasing wave steepness (vertical from top to bottom) the accuracy of the linear method decreases whereas the HOSM and NLSE₂ are also accurate for steeper sea states. For the steepest investigated sea states, the HOSM prediction is still on the accuracy level of the waveTUB input snapshot over very large distances, whereas the NLSE₂ shows to be less accurate. This indicates that the area of application of a forecast tool can be extended significantly by applying more complex wave models. In this context, the HOSM seems to be even applicable for very steep sea states.

Comparing the SSP of NLSE_FD (Figure 3) and linear simulation (Figure 4) shows no distinct trend. For the smallest steepness, both methods yield the same accuracy (independently from enhancement factor $\gamma$). With increasing wave steepness, the accuracy of both methods is still similar with a marginal tendency for better prediction with the NLSE_FD.

The influence of the enhancement factor on the prediction accuracy of the HOSM, NLSE₂ and linear simulations cannot be determined clearly. No trend is visible for the irregular sea sates with the smallest wave steepness. The discrepancy between linear transformation and HOSM simulation seems to increase slightly for increasing enhancement factor. Hereby, it remains unclear if this is caused by the different accuracies of the input snapshots or due the increase of the enhancement factor. The influence of the radar shadow on the accuracy of the prediction does not show a clear trend. For most of the cases, the influence of the radar shadow is negligible, otherwise the input snapshot without radar shadow provide a slightly improve as well as worsen of the accuracy.

Analysing the accuracy of the waveTUB input snapshot (black plus) reveals that the accuracy of the reproduction decreases with increasing wave steepness. In this context, the diagrams show that the HOSM prediction based on the waveTUB input snapshot remain almost on the same accuracy level (of the waveTUB input snapshot) for all wave steepness. This illustrates clearly that the accuracy of a real world wave prediction tool depends strongly on the accuracy of the detected input wave sequence. Strictly speaking, the wave prediction methods necessary for an accurate prediction are already available and not only the results of this investigation (cf. [10,13,14,15,26,30]) reveals that depending on the specific task very complex and numerically less efficient methods are not necessarily required for a successful application. However, the detection of the surrounding wave field by applying sophisticated wave inversion algorithm on X-band radar clutter is still challenging, but also the crucial factor for a successful application of a wave forecast tool. Hereby, the principle of measurements of X-band radars yield a main drawback as the radar cannot detect wave troughs behind (steep) wave crests. The main challenge hereby is, that most of the relevant applications regarding wave prediction and decision support are related to small wave steepness, where the intensity of the radar clutter is significant smaller compared to more steeper sea states. Less intensity correlates with reduced accuracy, particularly in the far field of the radar reducing the wave prediction horizon. Consequently, the fast linear method accurate enough for small wave steepness is significantly affected by the radar measuring accuracy. In addition, the vessel motions have to be measured accurately (for cruising vessels) in order to determine the position of the X-band radar as a critical prerequisite for the wave inversion. Generally, also the measuring accuracy (radar clutter as well as vessel motion) plays a major role for accurate wave prediction which can be divided in systematic and random error. The mentioned parameters relevant for an accurate wave prediction are only summarized without any detailed discussion as this is out of the scope of this paper.

4. Conclusions

This paper presents a numerical and experimental study on the applicability and limitations of use of intended prediction methods of different complexities: linear wave solution, the NLSE, NLSE_FD, NLSE₂, and the HOSM. The focus lies on the investigation of the prediction accuracy for varying parameters such as wave steepness, enhancement factor of the JONSWAP spectrum and wave propagation distance. For this purpose, irregular sea states with varying parameters are addressed.

The irregular sea state investigations focusses on JONSWAP spectra with varying wave steepness and enhancement factor. In addition, the influence of the propagation distance as well as the forecast horizon is discussed. The accuracy of the predictions are evaluated quantitatively by the SSP. The results show that the linear method is sufficient for the prediction of irregular sea states with small steepness with the same accuracy as the more complex methods. For increasing wave steepness, non-linear effects are more dominant resulting in a significant decreasing of the accuracy of the linear method. The NLSE results have shown throughout that this method is not applicable for deterministic wave prediction of arbitrary irregular sea states compared to the other methods. However, inclusion of the full dispersion operator resulted in a significant increased accuracy compared to the classical first-order NLSE. Both the NLSE₂ and HOSM results illustrate that deterministic wave prediction can be extended to steeper waves by applying more complex wave models. Hereby the HOSM shows a high accuracy also for the steepest investigated sea states over very large distances. In this context, the HOSM prediction based on the waveTUB input snapshot remains almost on the same input snapshot accuracy level for all wave steepness.

The accuracy of the waveTUB reproduction is throughout adequate but showed also that the accuracy decreases with increasing wave steepness of the irregular sea states. This illustrates clearly that the accuracy of a real world wave forecast tool depends strongly on the accuracy of the detected input wave sequence, i.e., assuming a radar-based wave detection, an accurate wave inversion algorithm used to extract the deterministic wave field from the radar clutter is the crucial factor for a successful application of a wave forecast tool.