## Abstract

Wave statistical properties and occurrence of extreme and rogue waves in crossing sea states are investigated. Compared to previous studies a more extensive set of crossing sea states are investigated, both with respect to spectral shape of the individual wave systems and with respect to the crossing angle and separation in peak frequency of the two wave systems. It is shown that, because of the effects described by Piterbarg, for a linear sea state the expected maximum crest elevation over a given surface area depends on the crossing angle so that the expected maximum crest elevation is largest when two wave systems propagate with a crossing angle close to 90°. It is further shown by nonlinear phase-resolving numerical simulations that nonlinear effects have an opposite effect, such that maximum sea surface kurtosis is expected for relatively large and small crossing angles, with a minimum around 90°, and that the expected maximum crest height is almost independent of the crossing angle. The numerical results are accompanied by analysis of the modulational instability of two crossing Stokes waves, which is studied using the Zakharov equation so that, different from previous studies, results are valid for arbitrary-bandwidth perturbations. It is shown that there is a positive correlation between the value of kurtosis in the numerical simulations and the maximum unstable growth rate of two crossing Stokes waves, even for realistic broadband crossing sea states.

## 1. Introduction

In many ocean regions sea states often consist of more than one wave system, so that the wave spectrum has multiple peaks with different directions of propagation and different peak wave periods. Such sea states are often referred to as crossing seas and are typical in the case that a local wind sea coexists with a system of swell.

From a sea captain’s point of view, a crossing sea may be a more challenging situation, since waves may hit the ship from different directions simultaneously. Crossing sea states have also attracted some attention in the context of rogue and extreme waves. Rogue waves (also called freak or abnormal waves) are waves that are unlikely to occur in a given sea state, given the averaged properties of a sea state, typically the significant wave-height *H*_{s}. In the last decades there has been a considerable interest in rogue waves, with a particular focus on understanding the mechanisms behind the formation of such waves as well as understanding in what types of sea conditions such waves may occur more frequently. Recent reviews on oceanic rogue waves can be found in Adcock and Taylor (2014) and, with a particular focus on applications to the marine industry, in Bitner-Gregersen and Gramstad (2016).

In the absence of inhomogeneities such as currents and variable bottom topography, the main mechanism that has been suggested to explain the occurrence of rogue waves beyond what is expected from standard linear or second-order wave theories is the effect of modulational instability. In its basic form the modulational instability is the phenomenon that a uniform wave train (a Stokes wave) is unstable to side-band perturbations, as first shown by Benjamin and Feir (1967). This instability, which may result in the formation of large individual waves, may also occur in the more realistic case of irregular wave fields described by a finite-width wave spectrum. It is, however, shown that the modulational instability is suppressed for spectra with broad bandwidth or wide directional spreading (Alber 1978; Crawford et al. 1980). This is in agreement with both numerical (Onorato et al. 2001; Gramstad and Trulsen 2007; Xiao et al. 2013) and experimental (Onorato et al. 2009, 2004) studies, which have shown that the occurrence probability of rogue waves is reduced when the wave spectrum becomes broader and has a wider directional spreading. Based on such results it has recently been an increasing concern whether effects of modulational instability are still relevant in realistic ocean conditions (Christou and Ewans 2014; Fedele et al. 2016; Benetazzo et al. 2015, 2017).

In this context, crossing seas have been suggested as a situation that may be associated with increased probability of rogue waves, even in realistic ocean conditions, and there are a number of previous studies addressing this topic, both in terms of modulational instability of systems of two crossing wave trains and in terms of numerical simulations of crossing random wave fields.

The modulational instability of a system of two Stokes waves with different directions of propagation has previously been considered within the framework of coupled nonlinear Schrödinger (NLS) equations. Roskes (1976) studied the modulational instability of a system of two uniform waves with the same frequency but different directions of propagation using a set of coupled NLS equations. Assuming perturbations along the mean propagation direction of the two Stokes waves, he showed that such a system is unstable for crossing angles 0° < *θ* < 70.5° and 136.1° < *θ* ≤ 180°. Based on this result it was suggested by Onorato et al. (2006) that more rogue waves may be expected in crossing sea states with crossing angles smaller than about 70°. Similar analysis was presented also in Onorato et al. (2010), who, based on the growth rate and amplification factor of an Akhmediev breather solution, suggested that a maximum amount of rogue waves could be expected for crossing angles between 20° and 60°. The more general case of two-dimensional perturbations was considered by Shukla et al. (2006), and stability analysis based on a more accurate fourth-order coupled NLS equation was derived in Gramstad and Trulsen (2011), who also considered the more general case where the two Stokes waves both have different directions of propagation and different frequencies. The effect of finite spectral width on the modulational instability of crossing systems was investigated by Shukla et al. (2007), who showed that finite bandwidth suppresses the modulational instability, just as in the case of unimodal seas (Alber 1978; Crawford et al. 1980).

The results from these existing studies indicate that a coupled system may experience increased instability growth rates compared to a single Stokes wave, in particular when the difference in propagation directions of the two waves is small. However, the instability results based on NLS-type equations suffer from nonphysical extensions of the instability regions beyond the bandwidth constraints of the equations, and it may be difficult to distinguish the physical instability regions with nonphysical finite-bandwidth effects (see, e.g., Gramstad and Trulsen 2011).

In this paper we derive the modulational instability of a system of two Stokes waves using the Zakharov equation (Zakharov 1968). The Zakharov equation has no bandwidth constraints, and the instability results are therefore valid for a wider range of configurations than existing results derived using NLS-type equations. The instability analysis of a single wave train using the Zakharov equation was first presented in Crawford et al. (1981) as well as in the following review by Yuen and Lake (1982). The modulational instability of standing waves (i.e., two counterpropagating waves with the same wavelength) was previously considered within the framework of the Zakharov equation by Okamura (1984). The derivation in this paper roughly follows the procedure of Crawford et al. (1981) and Okamura (1984) but with the necessary modifications to take the presence of two wave trains with general wavenumbers **k**_{a} and **k**_{b} into account.

Numerical simulations of two nearly unidirectional crossing wave systems, using a higher-order spectral method (HOSM), was presented in Onorato et al. (2010), who showed that the numerical simulations were in qualitative agreement with the results from the analysis of modulational instability, in the sense that more rogue waves were observed for crossing angles in the range of 40°–60°. Later, laboratory experiments of crossing unidirectional wave fields (Toffoli et al. 2011) also indicated an increased sea surface kurtosis for crossing angles of about 40°.

While the works mentioned above considered the crossing of two unidirectional or nearly unidirectional wave systems, numerical simulations of crossing sea states with more realistic directional spread wave spectra were performed by Bitner-Gregersen and Toffoli (2014) using HOSM. They showed that similar, although less pronounced, dependence on crossing angle could be identified also when the two wave systems were represented by more realistic spectra. Consistent with the theoretical and numerical results discussed above, it was shown by Cavaleri et al. (2012) that the accident of the cruise ship *Louis Majesty*, which was hit by a rogue wave that destroyed windows on deck 5 and caused two fatalities, which took place in a crossing sea state with crossing angle of about 40°–60°, and where the two wave systems had similar peak frequencies. A similar study on the sea state during the *Prestige* accident (Trulsen et al. 2015) showed that this accident also happened in a crossing sea state, with a crossing angle of about 90°. They were however not able to identify any increased probability of rogue waves for this sea state, also somewhat consistent with previous results mentioned above.

Note that the investigations that have indicated increased occurrence of rogue waves in crossing sea states have studied two wave systems with the same peak frequency and the same energy (same *H*_{s}). The opposite case that two crossing wave systems are well separated in frequency, typical in crossing seas consisting of wind sea and swell, was for example considered in Gramstad and Trulsen (2010), who in numerical simulations of an NLS-type equation, found no evidence of increased rogue wave occurrence due to the swell.

Since most of the studies on crossing seas have considered either the case that two wave systems have the same peak frequency and same energy or the case that two wave systems are very well separated in frequency, it is not clear to what extent existing results are robust also for more realistic cases, both with respect to spectral shape of the two crossing systems and with respect to the differences in peak frequency and energy of the two systems. The main objective of this paper is to investigate a somewhat wider range of crossing sea states, with different types of spectra and by letting two wave systems have different peak frequencies and different significant wave-height *H*_{s}. This is investigated both through theoretical analysis of the modulational instability of two crossing Stokes waves and by performing a large number of numerical simulations using the numerical wave model HOSM to simulate general broadband nonlinear sea states. In addition, linear effects on the statistical properties of crossing sea states are also considered.

The paper is organized as follows. First, in section 2 the analysis of the modulational instability of two crossing Stokes waves is presented. Then, in section 3 the more realistic case of random crossing sea states are considered, both in the context of linear theory, taking the effects described by Piterbarg (Piterbarg 1996; Krogstad et al. 2004) into account, and by numerical simulations of nonlinear phase-resolving numerical simulations. The relationship between the numerical results and the theoretical results from the instability analysis is discussed in section 4. Final discussion and conclusions are given in section 5.

## 2. Modulational instability of two crossing Stokes waves

As mentioned in the introduction, it has previously been suggested that there is a link between the enhancement in modulational instability experienced by a system of two waves and the increased occurrence of rogue waves in crossing sea states. These suggestions have been based on modulational instability results obtained from coupled NLS equations that suffer from bandwidth constraints. To obtain more accurate results of the modulational instability of two Stokes waves, we here derive the modulational instability of two crossing Stokes waves using the framework of the Zakharov equation. This represents a generalization of the modulational instability of a single Stokes wave, derived from the Zakharov equation by Crawford et al. (1981), and an extension of the bandwidth constraint results obtained from coupled nonlinear Schrödinger equations by Onorato et al. (2006), Shukla et al. (2006), and Gramstad and Trulsen (2011).

The analysis of modulational instability of two crossing Stokes waves using the Zakharov equation was previously carried out by Okamura (1984). Compared to the present analysis, the less general case of standing waves was considered; that is, **k**_{a} = −**k**_{b}, where **k**_{a} and **k**_{b} are the wavenumbers of the two Stokes waves. The analysis is a generalization of the results Okamura (1984) for general choices of **k**_{a} and **k**_{b}.

### a. Derivation of equations for instability

For a system of *N*_{m} discrete wave modes defined in terms of the complex amplitude spectrum *B*(**k**_{j}), *j* = 1, …, *N*_{m}, the Zakharov equation (Zakharov 1968) can be written in the form

Here, * denotes complex conjugate, ω is the angual frequency, and Δ*ω* = *ω*_{p} + *ω*_{i} − *ω*_{j} − *ω*_{m}. The indices *n, p, q*, and *r* run from 1 to *N _{m}*. Subscripts denote functional dependency of wavenumbers; for example,

*B*

_{j}=

*B*(

**k**

_{j}),

*T*

_{n,p,q,r}=

*T*(

**k**

_{n},

**k**

_{p},

**k**

_{q},

**k**

_{r}), and

*δ*

_{n,p,q,r}is the Kronecker delta

The exact expression for the Zakharov kernel function *T* can be found in, for example, Krasitskii (1994) or Janssen and Onorato (2007).

Extending the analysis of Crawford et al. (1981), we now consider a wave system consisting of two Stokes waves with wavenumbers **k**_{a} and **k**_{b} with corresponding complex amplitude functions *A*(*t*) and *B*(*t*), as well as a set of perturbations **k**_{a} ± **K** and **k**_{b} ± **K** with associated amplitudes *A*_{±}(*t*) and *B*_{±}(*t*), respectively. Inserting this into the Zakharov equation [(1)] and keeping only linear terms in the perturbation amplitudes gives the following system of six equations:

where subscripts *a*± denote **k**_{a} ± **K** and *b*± denote **k**_{b} ± **K**, and where the following short-hand notation is introduced, for example,

In the derivation of (3) we have made use of the symmetry properties of the Zakharov kernel function; that is, *T*_{a,b,c,d} = *T*_{b,a,c,d}, *T*_{a,b,c,d} = *T*_{b,a,d,c}, and *T*_{a,b,c,d} = *T*_{c,d,a,b}. Note that, as expected, the system of equations in (3) is invariant to the interchange of *A* and *B*.

One can show that (3) admits solutions in the following form:

where *Ω* is an angular frequency associated with the perturbations. Inserting (5) into (3) yields a set of four algebraic equations for the amplitudes *a*_{±} and *b*_{±}. One can show that this system of equations can be written as **x** = 0, where and where

where

For this system of equations to have nontrivial solutions the determinant must be equal to zero, det() = 0. This gives a fourth-order equation for Ω, and the corresponding unstable growth rate for the perturbations is given by ImΩ. We note that by setting *b* = *b*_{+} = *b*_{−} = 0 in the above equations the results of Crawford et al. (1981) for a single Stokes wave are recovered.

### b. Results of the instability analysis

In the following we consider two Stokes waves propagating with an angle *θ* relative to each other and with a possible difference in amplitude and wavenumber. To limit the number of free parameters, we restrict ourselves to the case that the two Stokes waves have the same wave steepness. That is, letting **k**_{a} and *a*_{a} and **k**_{b} and *a*_{b} be the wavenumber vectors and amplitudes of the two Stokes waves, we let

for *θ* ∈ [0, *π*] and *s* ∈ (0, 1]. For the sake of comparison with the numerical results presented in the next section, we let = 0.04 m^{−1} and = 1.414 m, where and are the peak period and significant wave height of wave system A chosen in the numerical simulations in the next section. The corresponding steepnesses of the two crossing Stokes waves are consequently *a*_{a}|**k**_{a}| = *a*_{b}|**k**_{b}| = 0.057. The effects of the wave steepness on the instability are discussed both in Crawford et al. (1981) for a single wave train and in Okamura (1984) for standing waves. Generally, increased steepness leads to enhanced growth rates and enlargement of the instability regions. However, the general features of the instability seem to remain, and in the following we consider only one value of the wave steepness.

Examples of instability regions in the **K**–plane are shown in Figs. 1–4. Generally, the instability regions of the combined system consist of regions that resemble the instability regions of the two waves separately (see also figures in Crawford et al. 1981), but in addition new regions of instability are also present, and in several cases these new regions have higher growth rates than the original instability regions. Note also that the instability regions are qualitatively similar to the ones derived from the coupled NLS-type equations in Shukla et al. (2006) and Gramstad and Trulsen (2011). Different from the coupled NLS results, however, results also outside the narrow-bandwidth regions can be trusted. For the case of standing waves (*θ* = 180° and *s* = 1), we reproduce the results of Okamura (1984). As shown by, for example, Fig. 3, the additional figure-eight-like instability region that was found by Okamura (1984) for standing waves is present also for not purely counterpropagating waves.

Figure 5 shows the normalized maximum growth rates as functions of the angle *θ* between the two waves for different values of the wavelength- and amplitude-ratio *s*. It is interesting to see that for relatively small angles there is a significant increase in the modulational instability compared to each of the wave systems alone. For angles close to 90°, however, there is almost no change in the instability compared to the unimodal case. This is consistent with previous studies on crossing sea states with crossing angle close to 90° (Trulsen et al. 2015) as well as with the numerical results in this study. For larger *θ*, a new figure-eight-like instability region with increased growth rates emerges. It is, however, worth mentioning that in this case, the growth rates of the “original” instability regions are actually slightly decreased, as shown by the dashed lines in Fig. 5.

It is also interesting to note from Fig. 5 that in the case the wavelengths of the two waves differ (*s* < 1), there is a weaker increase in the unstable growth rate. However, somewhat surprisingly, the enhanced modulational instability persists even for *s* = 0.5. For *s* = 0.25, however, the interaction between the two wave systems seems to be minor, and the growth rate is similar to that of one wave only, independent of the crossing angle. These results suggest that even typical wind-sea and swell systems may interact to increase the modulational instability, provided that the difference in wavelength of the two systems is not too large.

## 3. Statistical properties of random bimodal wave fields

While the analysis of modulational instability of two crossing Stokes waves presented above indicates that a system of two waves may experience enhanced unstable growth rates, it is not clear to what extent such theoretical results are relevant also for the more realistic case of random, short-crested bimodal wave fields in nature. In the following, linear and nonlinear numerical simulations of short-crested bimodal wave fields are performed, with the objective of investigating statistical properties and occurrence of rogue waves in crossing seas.

### a. Simulation setup

The numerical simulations in this study have been carried out using a numerical solver based on the HOSM, first proposed by Dommermuth and Yue (1987) and West et al. (1987). Short-crested bimodal wave fields have been simulated in a quadratic spatial domain with periodic boundary conditions. The nonlinear order *M* in the HOSM simulations was in this study set to *M* = 3, which includes the leading-order nonlinear dynamical effects, including the effect of modulational instability.

The horizontal plane is discretized using *n*_{x} × *n*_{y} = 512 × 512 grid points, and the maximum resolved wavenumbers were chosen as , where and are the peak wavenumbers of the two wave systems in the crossing sea state. Thus, the shortest wave that is resolved by the simulations, corresponding to the Nyquist frequency, spans two grid points, and the corresponding domain in the physical space covers an area of . For waves with *T*_{p} = 10 s on infinite water depth this corresponds to a computational domain of 5 km × 5 km. Note that these values are for the fully dealiased grid; the corresponding values before dealiasing when *M* = 3 are *n*_{x} × *n*_{y} = 1024 × 1024 and .

The wave fields are simulated in time for a total duration of *t*_{max} = 1800 s using a fixed-step ODE solver with integration time step Δ*t* = 0.3125 s, corresponding to Δ*t* = *T*_{p}/32 for *T*_{p} = 10 s. A weak dissipation of high wavenumbers is included to model the energy dissipation due to wave breaking, using the dissipation model suggested in Xiao et al. (2013). Some testing has been carried out to ensure that neither the choice of time step Δ*t* nor the dissipation model influence the statistical results significantly.

In all simulations the initial condition was chosen as a system of two JONSWAP spectra with cos^{N}(*ϕ* − *ϕ*_{p}) directional distribution with different peak directions of propagation and and different peak periods and . Random phases and amplitudes were assigned to the initial spectrum in all cases. Hence, each of the wave systems were chosen according to a spectrum of the type *E*(**k**) = *F*(*k*)*D*(*ϕ*), where **k** = (*k*_{x}, *k*_{y}) = *k*(cos*ϕ*, sin*ϕ*) and where

and

Here, Γ is the gamma function and the parameter *σ*_{A} has the standard values 0.07 for *k* ≤ *k*_{p} and 0.09 for *k* > *k*_{p}. The other spectral parameters *α*, *γ*, *k*_{p}, *ϕ*_{p}, and *N* were chosen to give the desired spectral shape, significant wave-height *H*_{s}, peak period *T*_{p}, and peak direction *ϕ*_{p}. Provided that the two wave systems are uncorrelated, the spectrum of the total wave field is *E*(**k**) = *E*_{a}(**k**) + *E*_{b}(**k**), where *E*_{a}(**k**) and *E*_{b}(**k**) are the individual spectra for the two crossing wave systems. The total significant wave height of the wave field can be defined as , while the significant wave heights of the individual systems are and . If the wave fields are uncorrelated, . Note that the energy-based definition of the significant wave height used in this paper (often denoted *H*_{m0}) differs from the definition based on the average of the 1/3 largest zero-crossing waves in a wave record (often denoted *H*_{1/3}). For linear and narrowband wave fields, these two definitions are almost equivalent. However, for broadband wave fields like crossing seas, *H*_{m0} and *H*_{1/3} will be different. However, in numerical simulations as carried out in this paper, *H*_{s} basically defines the energy of the initial spectrum, and an energy-based definition is most natural.

The spectral parameters for the different cases considered in this study are summarized in Table 1. For each case, different values of the angle *θ* between the two wave systems are considered, as well as different choices for the peak period and significant wave-height of the second wave system. However, as for the analysis of modulational instability in section 2b, in all cases and are chosen so that the mean wave-steepness is the same for both wave systems. Hence, for a given 0 < *s* ≤ 1, the period and significant wave heights are chosen so that . For the case of two identical wave systems (*s* = 1), 17 equally spaced values of *θ* between 0° and 180° are considered. For the cases of two wave systems with different peak periods and different significant wave heights (*s* = 0.75, 0.5, and 0.25), the five angles are considered. The peak directions and are chosen so that and , for .

The numerical simulations provide the full surface elevation *η*(*x*, *y*) at every integration time step (i.e., about every 0.3 s), and the solution was stored every 60 s, providing 31 snapshots of the surface in the spatial domain of 5 km × 5 km. All sea states were simulated with 20–50 repetitions with different initial random phases and amplitudes each time. Hence, in total this corresponds to an area of about 500–1250 km^{2} of wave data for each output time.

The focus of the present study is the occurrence of rogue waves, and in the following we discuss the following two statistical parameters that provide information about the occurrence of extreme and rogue waves in a wave field: the sea surface kurtosis *κ*_{4} and the maximum observed crest height *η*_{max}. That is, for a surface snapshot from the numerical simulations *η*_{i,j} = *η*(*x*_{i}, *x*_{j}) where *i* = 1, …, *n*_{x} and *j* = 1, …, *n*_{y} (here *n*_{x} = *n*_{y} = 512), we can define

where is the mean of *η*_{i,j}. Note that by construction of the numerical simulations, . In the presentation of results below, the values for *κ*_{4} and *η*_{max} are calculated as averages over all random realizations of the same sea state, and the maximum observed values over the duration of the simulations are plotted.

### b. Linear results

As a background for nonlinear numerical simulations presented in the following it is useful to briefly discuss the statistical parameters *κ*_{4} and *η*_{max} in the case of a purely linear random sea surface, in which case the surface is Gaussian distributed and can be viewed as a Gaussian random field.

It is well known that the kurtosis of a Gaussian sea surface is equal to three. Deviation from three is used as a measure of the degree of non-Gaussianity of the sea surface, in the sense that a kurtosis larger than three indicates more extreme waves than expected according to linear wave theory. However, one should keep in mind that *κ*_{4}, as defined in (10), is not the “real” population kurtosis, but the sample kurtosis of a realization of a Gaussian random field, and the sample kurtosis may be a biased estimator for the real population kurtosis. In fact, even for independent and identically distributed (IID) normal samples the expected value of the standard estimator for the sample kurtosis is biased so that *E*[*κ*_{4}] = 3(*n* −1)/(*n* + 1), where *n* is the number of samples (here *n* = *n*_{x}*n*_{y} = 262 144). However, as further discussed in the appendix, samples from a spatial snapshot of the surface elevation are not IID but contain values that are statistically dependent according to the relevant cross correlation of the sea surface, which in turn depends on the wave spectrum through the Wiener–Khinchin theorem. This dependency has the effect that the bias of the kurtosis estimator is significantly enhanced (Bai and Ng 2005), and deviation from the theoretical value of three must be expected, even for a Gaussian field. Note that the alternative formula for sample kurtosis that is unbiased for IID samples also suffers from this bias. Moreover, the magnitude of the bias depends on the cross correlation of the field so that a more narrow wave spectrum will give a stronger bias. However, as shown in the appendix, this effect is most pronounced for small simulation domains, so for the present numerical results, which applies to a relatively large computational area, this effect is of limited importance, but it highlights the importance of having a sufficiently large computational domain.

The theoretical treatment of the maximum wave elevation *η*_{max} is more complicated, even in the linear case. When considering a finite surface area, as in the present numerical simulations, the expected value for *η*_{max} naturally depends on the characteristics of the sea state, in particular on the typical number of waves in the domain. This, in turn, depends on the average values for the wave and crest length of the wave field. These effects are discussed in the context of Gaussian random fields in general by Piterbarg (1996), and in the context of ocean surface waves by, for example, Krogstad et al. (2004). It can be shown that the typical number of waves *N*_{w} in a rectangular domain of size *V* can be estimated by

where and are the mean wave and crest lengths, respectively, which can be derived from the wave vector spectrum *E*(**k**) as explained in Krogstad et al. (2004). It can further be shown that the expected value for *η*_{max} can be estimated as

Note that the Piterbarg theory, on which the above expressions are based, considers the asymptotic extreme value distributions of homogeneous Gaussian fields. Hence, the results based on the Piterbarg theory must be viewed as an approximation, which is approximately valid for linear wave fields.

For crossing sea states, it is clear that and depend on the crossing angle of the two wave systems, which will affect the expected value for *η*_{max} through the relation in (12). In particular, it is obvious that the typical crest length in the wave field will be shorter for crossing angles close to 90°. This effect is shown in Fig. 6, which shows the expected value *E*[*η*_{max}/*H*_{s}], according to the Piterbarg theory, for the different types of crossing sea states considered in this study. Thus, from purely linear considerations a higher maximum wave elevation is expected for crossing sea states with crossing angles close to 90°, while the lowest values for *η*_{max} are expected for co- and counterpropagating wave systems. It is also clear from the figure that the expected value of *η*_{max}/*H*_{s} is reduced when the two wave systems differ in peak frequency and energy (*s* < 1). The theoretical results are qualitatively confirmed by linear numerical simulations, shown by solid lines with symbols in Fig. 6. The quantitative difference between the Piterbarg theory and the linear simulations is due to the fact that the Piterbarg theory is, as mentioned above, based on approximations as well as due to sampling variability in the numerical results obtained from a limited number of simulations.

### c. Results from nonlinear simulations

Keeping in mind the results from the theoretical instability analysis of section 2, as well as the linear results of the previous section, we now turn to the nonlinear case and present results from the numerical simulations using the HOSM model. The main focus is how *κ*_{4} and *η*_{max} depend on both the angle *θ* between two crossing wave systems and the ratio of the wavelengths and significant wave heights of the two systems, as well as on other properties of crossing sea states.

To visualize how the statistical properties depend simultaneously on *θ* and *s*, Gaussian process (GP) regression has been used to fit surfaces in the *θ*–*s* plane representing the best fit to the simulated values of *κ*_{4} and *η*_{max}. This approach is useful both with respect to visualization of the results and also for estimating the uncertainty related to lack of observations in regions of the *θ*–*s* plane where few simulations are available. The GP regression is performed so that the resulting fitted surfaces goes exactly through the actual values obtained from the numerical simulations. Similarly, corresponding surfaces are fitted to the 95% confidence intervals, which are calculated using the standard deviation over runs with different random seeds.

Figure 7 shows how *κ*_{4} and *η*_{max} depend on the angle between the two wave systems, for the case that the two systems have the same peak frequency (i.e., *s* = 1), for the three different types of crossing sea states listed in Table 1.

It is evident from Fig. 7a that the kurtosis is larger for relatively small and relatively large angles and attains a minimum around *θ* = 90°. Note also that, consistent with previous works on crossing seas (Bitner-Gregersen and Toffoli 2014) as well as well-known results for unimodal wave systems (e.g., Gramstad and Trulsen 2007; Janssen 2003; Mori and Janssen 2006; Mori et al. 2011), kurtosis is higher for the most narrow spectrum (*γ* = 6, *N* = 100) than for the more broad and directionally spread cases. One should also notice the local peak in kurtosis at approximately *θ* = 45° for the case *γ* = 6, *N* = 100. This is remarkably consistent with previous works (Onorato et al. 2010; Toffoli et al. 2011) that have suggested that a maximum probability of rogue waves is obtained for 40° < *θ* < 60°. It was suggested by Onorato et al. (2010) that this feature could be explained from the modulational instability of crossing waves and that this range of angles represents a compromise between large growth rates and large amplification factors of the instability. Similar results have also been obtained previously from numerical simulations of two nearly unidirectional crossing wave systems. The present results indicate that this feature is relevant only for relatively narrowbanded crossing systems. However, one should note that it is likely that this conclusion may also depend on the wave steepness of the two crossing wave systems, which is not investigated here.

For the maximum crest height *η*_{max} shown in Fig. 7b, no clear dependence of the crossing angle or differences between the different types of spectra are observed. When interpreting this result in light of the linear results discussed above, however, it is clear that the nonlinear contribution has an opposite effect compared to the linear effects shown in Fig. 6. This is confirmed by Fig. 7c, which shows , where is the expected value of *η*_{max} according to the linear Piterbarg theory [(12)], shown in Fig. 6. Thus, nonlinearities seem to have a similar effect on *η*_{max} as observed for the kurtosis, but since the linear results due to Piterbarg theory are to a large extent opposite to the effects of nonlinearity, the actual observed *η*_{max} in the nonlinear simulations is almost independent of the crossing angle.

For all three types of crossing seas shown in Table 1, several cases were run for which the two crossing systems have different peak frequency and significant wave height in addition to different propagation angles. The results from these simulations are shown in Fig. 8, which shows how *κ*_{4} and *η*_{max} depend on the wavenumber and significant wave-height separation ratio for different crossing angles.

The general feature of Fig. 8 is that for the two most narrow spectra (cases B and C in Table 1) both the kurtosis and *η*_{max} show a decreasing trend when the two wave systems become more separated in frequency space (corresponding to decreasing values for *s*). However, for several crossing angles the magnitude of both kurtosis and *η*_{max} show almost no change for *s* ≥ 0.5. For the largest separation (*s* = 0.25), however, lower values for kurtosis and *η*_{max} are consistent for most angles. This indicates that the nonlinear effects shown in Fig. 7, which tend to increase the occurrence of rogue waves, is relevant also for wave systems that are moderately separated in frequency space. Note also that, as shown in Fig. 6b, both the effect of nonlinearity and the linear effects due to Piterbarg theory have the same consequence, namely, that increasing the separation frequency tends to decrease the observed value of *η*_{max}.

For the most broadbanded case (*γ* = 1, *N* = 4), the picture is less clear, and both the kurtosis and show a flat, and in some cases even increasing, trend for decreasing values of *s*. However, *η*_{max}/*H*_{s} is still decreasing with decreasing *s*, likely because of the fact that the linear effects dominate in this case.

## 4. Comparison of instability analysis and results from numerical simulations

Given the results of the previous sections, it is now of interest to discuss the statistical results from the numerical simulations in light of the theoretical stability analysis in section 2. In fact, it is interesting to observe that the curves of unstable growth rate in Fig. 5 share some common characteristics with Figs. 7a and 7c, showing the kurtosis and normalized maximum crest elevation , in the sense that these quantities are largest for small and large separation angles and have a minimum for angles around 90°. Although this similarity seems to be most pronounced for the most narrow spectrum (*γ* = 6, *N* = 100), some resemblance is seen also for the two broader cases.

This is confirmed in Fig. 9, which shows the relation between the unstable growth rate from the instability analysis and the statistical results from the numerical simulations, where the growth rate is taken from the case of two Stokes waves with the same crossing angle and the same amplitude ratio as used in the corresponding numerical simulations (i.e., same values for *θ* and *s*). It is clear from the figure that for the two most narrow cases (cases B and C in Table 1) there is a positive correlation between the unstable growth rate and both the kurtosis and the maximum surface elevation normalized by the expected value from the linear theory . However, because of the effects discussed in section 3b, no clear correlation is observed between the maximum growth rate and the actual observed maximum crest height *η*_{max}/*H*_{s}. These conclusions are also confirmed by Fig. 10, which shows the corresponding correlation coefficients together with the *p* value for the statistical test that the correlation is different from zero.

For the most broadband case (*γ* = 1, *N* = 4) there is no statistically significant correlation between the growth rate and the kurtosis or , indicating that, as one might expect, as the spectrum becomes broader, the effect of the modulational instability is suppressed. This feature is well known in unimodal sea states, and it is not surprising that this also applies to crossing seas, which is also consistent with the theoretical work of Shukla et al. (2007).

One should also note that for the two broadest cases (cases A and B in Table 1), there is a positive correlation between the unstable growth rate and *η*_{max}/*H*_{s}. This is due to the fact that the Piterbarg theory predicts decreasing *η*_{max} with decreasing *s* (see Fig. 6), and similarly, the modulational instability predicts decreasing growth rate for decreasing *s* (see Fig. 5). However, since the Piterbarg effect is a purely linear effect while the modulational instability is a nonlinear effect, these two effects are not related, and this correlation does not represent a causal relationship.

## 5. Conclusions

The motivation for this work is several previous investigations suggesting that crossing sea states may lead to increased probability of rogue waves. In particular, it has been suggested that more rogue waves may be expected if two wave systems with similar peak frequency and similar energy propagate with crossing angle close to 50°.

In this paper, several types of crossing sea states have been simulated numerically, investigating the effects of crossing angle, differences in peak frequency and energy, and the influence of the spectral shape of two crossing systems. For narrowband and almost unidirectional crossing wave systems with the same peak frequency and the same energy, our numerical results confirm an increased sea surface kurtosis for crossing angles close to 50°, although even higher values of kurtosis are observed for very small crossing angles, where the resulting total spectrum essentially is a unimodal narrow spectrum. Minimum kurtosis is observed for crossing angles close to 90°, with larger values for both smaller and larger separation angles.

Similar results are also obtained for more realistic broadband and directionally spread crossing wave systems, which also show a minimum kurtosis for crossing angles close to 90° and largest values for small and large crossing angles. However, in these more-broadband cases, the local peak in kurtosis for crossing angles of about 50° does not appear.

Further, the nonlinear simulations show that kurtosis is generally decreasing when two wave systems become more separated in frequency space. However, for relatively weak separation in frequency space the kurtosis is only slightly reduced.

Interestingly, the dependence of kurtosis on crossing angle and frequency difference in the numerical simulations is in qualitative agreement with the results from the analysis of the modulational instability of two crossing Stokes waves, in the sense that the largest growth rate is observed for large and small crossing angles, with a minimum around 90°, and that increasing difference in the frequency of the two waves results in smaller growth rates. This correlation between the maximum growth rate of two interacting Stokes waves and the simulated kurtosis is observed also for the relatively broadband crossing system with *γ* = 3.3, *N* = 16, but not for the broadest spectrum with *γ* = 1, *N* = 4. This suggests that for broader spectra the effect of the modulational instability is weakened; however, the results suggest that nonlinear modulational-instability effects due to the interaction of the two wave systems play a role also in relatively broadband, realistic types of crossing sea states.

For wave-statistical properties that depend on the average size of the individual waves in a wave field, such as the maximum crest height *η*_{max} in a given area, there is another effect that also plays a role in crossing sea states. Namely, the fact that the typical crest and wave lengths in a crossing sea state strongly depend on the crossing angle and frequency difference of the two wave systems, and hence influence the expected value of, for example, *η*_{max} even in the linear case. Interestingly, this effect acts in opposite direction to the nonlinear effects, so that *η*_{max} in a nonlinear crossing sea state seems to be almost independent of both the crossing angle and the spectral bandwidth of the crossing systems.

It should be mentioned that the nonlinear effects in crossing sea states depend on the average wave steepnesses of the two wave systems. In this paper only one value for the wave steepness of the two crossing wave systems is considered, with the restriction that the wave steepnesses of the two systems are equal. Hence, for an even more complete picture of the effect of crossing sea states on wave statistical properties, a larger set of steepnesses should be investigated, including the effect of varying the steepnesses of the two wave systems independently.

## Acknowledgments

This work is funded by the Research Council of Norway Grant 256466.

### APPENDIX

#### Effect of Domain Size on Estimation of Sea Surface Kurtosis

An effect that, to our knowledge, has to a large extent been overlooked in studies on statistical properties of waves is the fact that a realization of a sea surface in a limited spatial area (or a realization of a time series of a limited duration) contains a large number of dependent statistical samples compared to overall number of samples. In such a case, the expected values of the estimates of, for example, the kurtosis and skewness differ from their theoretical values for a linear sea surface. In fact, the effect of statistical dependence in a realization of a linear sea surface will tend to underestimate the kurtosis (Bai and Ng 2005; Bao 2013) because the estimators for kurtosis and skewness are biased for dependent data of limited sample sizes.

Figure A1a shows the estimated kurtosis of linear realizations of a sea state with peak period *T*_{p} = 10 s, calculated as a spatial average over two different domain sizes: a 5 km × 5 km area discretized with *n*_{x} × *n*_{y} = 512 × 512 points and a 1.25 km × 1.25 km area discretized with *n*_{x} × *n*_{y} = 128 × 128 points. For neither of the two domain sizes does the kurtosis converge exactly to the true value of three with increasing number of realizations, and for the smallest domain size, the kurtosis is underestimated by about 3%–4%. This is also illustrated by the plot in Fig. A1b, which shows that the distributions of the sample kurtosis are skewed with mean below three.

The same effect is also relevant for kurtosis calculated from time series of linear surface elevation, as illustrated in Fig. A1c, which shows the averaged kurtosis obtained from 100 000 random time series realizations of two different sea states. Consistent with the results in Bai and Ng (2005), the kurtosis is underestimated for time series of short duration and only slowly approaches three with increasing duration. Note also that, as shown in Bai and Ng (2005), the bias of the sample kurtosis will naturally depend on the serial dependence in the wave field. Consequently, a narrowband wave field with a long correlation length will give a stronger underestimation of the kurtosis compared to a more broadband wave field with a shorter correlation length. This effect is shown in Fig. A1c by the fact that the JONSWAP spectrum with *γ* = 7 has lower kurtosis than the more broadband spectrum with *γ* = 1.

Hence, one should be aware that the kurtosis, which is often used as an indicator for effect of nonlinearity in a sea state, may be underestimated in numerical simulation or field measurements when using a small computational domain or a short time series, purely because of the effect of bias in the estimators for kurtosis for dependent samples.

## REFERENCES

*Proc. 23rd Int. Conf. on Offshore Mechanics and Arctic Engineering*, Vancouver, BC, Canada, American Society of Mechanical Engineers, OMAE2004-51336, https://doi.org/10.1115/OMAE2004-51336.

*Asymptotic Methods in the Theory of Gaussian Processes and Fields*. American Mathematical Society, 206 pp.

## Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).