Wave Potential

Stokes’ theory of ocean waves assumes irrotational flow of an incompressible fluid, described by a scalar velocity potential \(\mathrm {\Phi(\vec{x},t)}\) which satisfies

with a velocity field \(\mathrm {\vec{u}}\) given by

Assuming an impermeable sea-floor of constant mean depth \(\mathrm {h}\), we have \(\mathrm {u_z(z=-h)=0}\), which we write in terms of the wave potential as

The free surface of the wave is given by \(\mathrm {z=\eta (x,y,t)}\). For sea states with long-crested waves the \(\mathrm {y}\) coordinate is irrelevant, while for short-crested waves it is sometimes useful to transform between a cartesian and a polar coordinate system, such that \(\mathrm {z=\eta (r,\theta ,t)}\). Now, \(\mathrm {r}\) plays the role which \(\mathrm {x}\) plays in long-crested waves, and \(\mathrm {\eta }\) results from the superposition of several \(\eta _{\mathrm i}(r,t)\equiv\eta (r,\theta _{\mathrm i},t)\).

Two boundary conditions are imposed at the wave surface. The kinematic boundary condition

which states that a particle at the surface moves together with the surface, and the dynamic boundary condition

which follows from assuming constant pressure along the surface. Here, \(\mathrm {g}\) is the acceleration of gravity.

Together with Equation (24), which applies within the bulk of the ocean, Equation (26), specify the set of differential equations from which \(\mathrm {\Phi}\) is obtained.

In order to obtain analytic results for \(\mathrm {\Phi}\), a perturbative approach is taken.

The velocity potential at the wave surface is extrapolated onto the mean surface level \(\mathrm {z=0}\) using a Taylor expansion, with terms of higher order in the expansion expected to contribute at higher order in the ordering parameter \(\mathrm {\epsilon}\).

The expansion is, to order \(\mathrm {\mathcal{O}(\epsilon^2)}\),

Likewise, the surface elevation is written as a series expansion,

where it is usually assumed that \(\mathrm {\eta ^{(1)}}\) is a linear wave. This is a construction which allows us to study the non-linear effects which appear at second order and beyond.

The surface elevation is a function of \(\mathrm {x}\), \(\mathrm {y}\) and \(\mathrm {t}\), or equivalently and for many purposes more conveniently, a function of \(\mathrm {r}\), \(\mathrm {\theta }\) and \(\mathrm {t}\). Staying with the former set of coordinates, the surface elevation in Fourier space is a function of wavenumber \(\mathrm {\vec{k}=(k_x,k_y)}\) and angular frequency \(\mathrm {\omega }\). They are related through the dispersion relation for surface waves at finite water-depth,

which is complete to second order. At third order, the dispersion relation includes an additional term.

The first order contribution \(\mathrm {\eta ^{(1)}}\) is provided to WaveKin by its Fourier components \(a_{\mathrm i}\),

where the sums run over the number of positive frequencies \(\mathrm {N}\) and the number of directions \(\mathrm {M}\) in the discretized directional spectrum, \(\psi_{\mathrm {ik}}=\omega _{\mathrm i}t-\vec{k_{\mathrm {ik}}}\cdot \vec{x}-\varepsilon _{\mathrm {ik}}\) is the phase function, with angular frequency $_{i}={i}$. The phase is \(\varepsilon _{\mathrm {ik}}\) and the wavenumber \(\vec{k_{\mathrm {ik}}}\) is obtained from \(\omega _{\mathrm i}\) through use of the dispersion relation and wave direction \(\theta _{\mathrm {ik}}\). The two subscripts denote a discretization of the Fourier components of a directional spectrum. In the following, the second subscript \(\mathrm {k}\) is dropped in favor of a compact notation. It should be implicitly understood that for directional spectra, all sums over Fourier components also imply a sum over directions.

For a simulated wave from a given spectrum, the random phase will typically be drawn from a uniform probability distribution. For a directional spectrum this involves \(\mathrm {MN}\) independent and identically distributed random numbers. See Section 4.2 for a further discussion on this.

The sum of Equation (1) can be re-written to run over positive frequencies only. Because \(\mathrm {\eta (x,y,t)}\) is real-valued, \(a_{-\mathrm i}=\bar{a}_{\mathrm i}\). This gives

where the sum over directions is now implicit.

Velocities and accelerations are also real-valued; WaveKin therefore works with one-sided spectra only. Time-series may have a non-zero mean value, therefore the one-sided spectra include the zero-frequency components. The zero-frequency component for surface elevation \(\mathrm {a_0}\) is also included in the one-sided spectrum, though its value should be zero or near zero. See Section 4.1.6 for how symmetries allow summations of first and second order terms to run over positive frequencies only and Section 4.1.7 for a further discussion on the zero-frequency components.

The expansion is in terms of an ordering parameter \(\mathrm {\epsilon}\). It can be seen that higher orders in the expansions of the boundary conditions at the wave surface results in terms with higher powers of spatial derivatives, surface elevation, and the velocity potential. In Fourier space, spatial derivatives result in additional factors of \(\mathrm {k}\), and we have constructed the first order surface elevation to be a sum of harmonics with amplitude \(\mathrm {|a|}\).

The combination of factors of \(\mathrm {k}\) and \(\mathrm {a}\) provides a physical meaning to the ordering parameter. The wave steepness is proportional to \(\mathrm {ak}\), for a harmonic wave elevation \(\mathrm {\eta _k=a\sin kx}\). In order for the series expansion to converge, the linear wave should have a spectral content which satisfies \(|a_{\mathrm i}|k_{\mathrm i}<\epsilon\) for all \(\mathrm {i}\).

In addition to the Taylor series expansion of the wave surface boundary conditions and the expansion of \(\mathrm {\eta }\), the velocity potential is written as a series expansion

A solution is obtained for each term of the expansion by ordering the terms in the wave surface boundary conditions, setting \(\mathrm {z=\eta }\) in Equation (1). Terms at equal order in \(\mathrm {\epsilon}\) are solved for independently.

Solving for the first order velocity potential gives

Following Marthinsen and Winterstein (1992) for long-crested waves and Marthinsen and Winterstein (1992) for short-crested waves, the velocity potential and second order contribution to surface elevation is given to second order by a set of rather complex expressions, organized in terms of sum-frequency and difference-frequency contributions. The expressions are analytic, and therefore support differentiation, allowing the analytic calculation of velocities and accelerations.

Solving for the second order surface elevation gives

Solving for the second order velocity potential gives

where the latter term is independent of spatial coordinates and therefore irrelevant for kinematics.

The second order kernels \(E_{\mathrm {i,j}}\) and \(P_{\mathrm {i,j}}\) are given by

and

The solutions for both \(\mathrm {\eta ^{(2)}}\) and \(\mathrm {\Phi^{(2)}}\) depend on products, sums and differences of Fourier components. Because the two-sided spectrum can be expressed as a one-sided spectrum, the second order contributions are naturally separated into sum-frequency and difference-frequency contributions.

Because the sum and difference-frequency components only depend on the sum and difference of frequencies (wavenumbers) and the Fourier components satisfy \(a_{-\mathrm i}=\bar{a}_{\mathrm i}\), the calculations of second order contributions to surface elevation, velocities and accelerations can be considerably simplified.

The kernels from Equation (39) and Equation (40) can be written in terms of sums and differences of frequencies and corresponding wavenumbers,

and

where \(\mathrm {i,j}>0\).

The symmetry properties of the second order kernels are

where \(\mathrm {i,j}>0\).

For the second order wave elevation of Equation (1) we now have \(\mathrm {\eta ^{(2)}=\eta ^{}\eta ^{-}}\) with

The second order velocity potential of Equation (38) can now be written as \(\mathrm {\Phi^{(2)}=\Phi^{}\Phi^{-}}\) with

where space-independent terms have been omitted.

The difference-frequency kernels \(E_{\mathrm {ij}}^-\) and \(P_{\mathrm {ij}}^-\) are singular for \(\mathrm {i}=\mathrm {j}\). If one chooses to take the analytic limit, one gets a non-zero contribution to the zero-frequency contribution of the difference-frequency transfer functions. For \(E_{\mathrm {ii}}^-\) this gives

where

is the group velocity.

Equation (49) follows Marthinsen and Winterstein (1992).

With non-zero contributions from the difference-frequency kernels when \(\mathrm {i}=\mathrm {j}\), the resulting transfer functions give finite contributions to the zero-frequency component of the second order elevation and velocity potential. This corresponds to finite mean-values for \(\mathrm {\eta }\), \(\mathrm {u_x}\) and \(\mathrm {u_y}\). This can optionally be avoided by setting all zero-frequency components to zero, (equivalently, setting the diagonal entries of the difference-frequency kernels to zero) however this procedure lacks physical or theoretical justification.

In particular, the so-called mean sea state set-down leads to an inconsistency with the boundary condition applied at the sea floor, i.e. at \(\mathrm {z=-h}\), as the sea floor depth \(\mathrm {h}\) should properly be measured from the mean sea level, \(\mathrm {<\eta >}\), not from the first order contribution to the same, \(\mathrm {<\eta ^{(1)}>}\). The inconsistency might be small for most practical uses, but deserves mentioning. This is discussed occasionally in the literature, e.g., Marthinsen and Winterstein (1992) and The default behaviour of WaveKin in RIFLEX is to set the diagonal terms equal to zero.

Depending on how the target spectrum is filtered to obtain the first order wave elevation defined by \(\{a_{\mathrm {i}}\}\), the sum-frequency transfer functions may result in significant high frequency contributions to the second order wave elevation and kinematics. The acceleration terms are particularly sensitive to the high frequency part of the spectrum. Care must be taken when defining the first order wave elevation to avoid un-physical high frequency terms arising at second order. Otherwise, un-physically large accelerations could result.

The issue of large accelerations, i.e., rapidly varying velocities, is most acute near the mean sea surface and above (\(\mathrm {z\ge0}\)).

By target spectrum is meant the empirically established power spectrum \(\mathrm {S(\omega )}\) for long-crested waves or the directional spectrum \(\mathrm {S(\omega ,\theta )}\) for directional sea states.

In order to generate an irregular sea state, phases are typically drawn from a uniform distribution, with the amplitude of Fourier components fixed by a filtered power spectrum. The resulting wave elevation is then referred to as a linear wave and identified as the first order contribution in the perturbative expansion of \(\mathrm {\eta }\).

Various heuristics are presented in the literature for how the filtering should be done, see e.g. Marthinsen and Winterstein (1992) and Marthinsen and Winterstein (1992), but no solid theoretical foundation underlies these methods. For many purposes, results may be insensitive to the details of the filtering, however results from WaveKin indicate that filtering has a significant effect on wave kinematics near the sea surface at water depths of \(20-40{\mathrm {~m}}\) and significant wave heights in the range of \(5-6{\mathrm {~m}}\).

Equating the linear wave with the first order wave elevation is a matter of convenience, but could under-estimate the non-linear nature of Stokes’ waves. The assumption of un-correlated phases is justified by very weak wave-wave interactions, which in turn applies to low wave steepness. While the first order wave elevation will always have a more narrow banded spectrum than the full wave, it will nevertheless contain steep waves.

A directional spreading function \(\mathrm {D(\theta ,\omega )}\) is defined by

together with the normalization condition

The power spectrum \(\mathrm {S(\omega )}\) for directional sea states is obtained by integrating over the directional coordinate

The simplest approach to describing a directional sea state is to assume a separable directional spectrum,

where the power spectrum \(\mathrm {S(\omega )}\) (e.g., a JONSWAP spectrum) and the directional spreading function \(\mathrm {D(\theta )}\) are given separately and together define the sea state.

A common choice for \(\mathrm {D}\) is a frequency-independent cosine-2s distribution

for \(\mathrm {|\theta -\theta _0|\leq\frac{\pi }{2}}\) and zero otherwise. \(\mathrm {C_s}\) is a normalization factor. The width of the distribution is parametrized by \(\mathrm {s}\) and the mean direction is given by \(\mathrm {\theta _0}\).

Frequently, convective terms are neglected in the calculation of wave kinematics. In the current implementation of WaveKin the inclusion of these terms is optional. If they are omitted, only the contribution from partial time-derivatives are included

Convection terms can optionally be included in the accelerations returned by WaveKin. These are calculated by taking spatial derivatives of the velocity field (in the frequency domain) and multiplying with the velocity field (in the time domain).

Convection terms contribute to the accelerations of momentum-carrying fluid particles, and thus contribute to the inertia term of Morison’s load model. They emerge when we take the total derivative of velocity with respect to time,

where \(\mathrm {\vec{v}}\) is the velocity vector and \(\mathrm {\vec{u}=\vec{u}(\vec{x},t)}\) is the velocity field.

The inclusion of convection terms is turned off by default in WaveKin. This follows an assumption that they are small and can be neglected for practical purposes. This assumption should be questioned for severe sea states. The last term in Equation (59) constitutes the convection terms. It should be seen that its magnitude depends on the gradient of the velocity field, meaning that rapidly varying velocities, as found near the sea surface in severe sea states, may lead to significant contributions from convection terms.

The velocity potential \(\mathrm {\Phi(x,y,z,t)}\) is extrapolated above \(\mathrm {z=0}\) according to the Taylor expansion of Equation (1). Kinematics are obtained from the extrapolated potential by spatial and temporal differentiation in the same manner as for the velocity potential at and below \(\mathrm {z=0}\). However, it should be noted that the second order contribution at \(\mathrm {z>0}\) only involves the first order velocity potential \(\mathrm {\Phi^{(1)}}\). There are no quadratic transfer functions involved above mean surface level until third order contributions.

High frequencies resulting from providing an improperly filtered spectrum as the input first order wave may cause unphysical accelerations near to and above the mean surface level. To compensate for this, an option is provided to filter out large frequencies prior to extrapolating the kinematics above the mean surface level.

Together with, or instead of this option, there is an option to set a threshold surface level \(z_{\mathrm {lim}}<0\), above which kinematics are extrapolated from the potential at \(z=z_{\mathrm {lim}}\). As the large frequencies are rapidly dampened even at very small values of \(\mathrm {z<0}\), setting \(z_{\mathrm {lim}}\) to only a few centimeters below \(\mathrm {z=0}\) may be sufficient to avoid the excessive spikes in acceleration.

The default behaviour of WaveKin in RIFLEX is to use the entire provided spectrum when extrapolating kinematics above the mean surface level, and to extrapolate from \(\mathrm {z=0}\).

In order to extrapolate kinematics between two nodes, one of which is dry, it can be impractical to set kinematics to zero above the sea surface \(\mathrm {z>\eta }\) (sometimes referred to as a dry node). Instead, the option exists to set kinematics at \(\mathrm {z>\eta }\) equal to the the kinematics at the surface: \(\mathrm {z=\eta }\). This is the default behaviour of WaveKin in RIFLEX.

For kinematics on the surface and at or above the mean surface level, \(\mathrm {z=\eta \ge0}\), the kinematics are obtained at second order by extrapolation of the first order kinematics at \(\mathrm {z=0}\). For kinematics on the surface and below the mean surface level, \(\mathrm {z=\eta <0}\), there is a choice between using the second order potential, which exists at all \(\mathrm {z\le0}\) or using the first order potential at \(\mathrm {z=\eta <0}\) together with the first order potential extrapolated from \(\mathrm {z=0}\). The latter method provides a consistent calculation of kinematics on the surface regardless of whether it is above or below \(\mathrm {z=0}\). However, it causes a discontinuity in the kinematics on and directly beneath the surface.

The magnitude of the discrepancy between the two methods should be of third order, if the series expansion behaves consistently. Thus, a comparison of the two methods is a test of the series expansion.