3.1. Assumptions

Considering the wave field as random, stationary and Gaussian, the wave prediction problem can be formulated in a rigorous statistical framework, as shown by Simanesew et al. (Reference Simanesew, Trulsen, Krogstad and Nieto Borge2017) and Mérigaud & Ringwood (Reference Mérigaud and Ringwood2018a ). The assumptions underlying the approach are commonplace in oceanography and marine engineering since the 1950s (Pierson Reference Pierson1952).

• (i) Assumption 1. The FSE $\eta$ is a zero-mean and Gaussian random field, which means that any finite set of wave elevation samples $\eta (\textbf{x}_1,t_1), \eta (\textbf{x}_2,t_2), \cdots \eta (\textbf{x}_N,t_N)$ , taken at $N$ positions and instants, follows a multivariate Gaussian distribution. The Gaussian nature of ocean waves is verified, at least to a good approximation, in the majority of sea conditions of interest for marine operations (Ochi Reference Ochi2005), which excludes major storms and shallow water waves (with marked peak-to-trough asymmetry), as well as specific phenomena such as wave breaking, whitecapping and rogue waves.

• (ii) Assumption 2a. Examined over the time scales of interest (typically, a few seconds to 10 min), the wave field can be considered stationary, i.e. its statistical properties do not change substantially. Since the waves are also Gaussian, the field stationarity reduces to the stationarity of its first- and second-order moments (i.e. the wave process is wide-sense stationary).

• (iii) Assumption 2b. The wave field is statistically homogeneous over the spatial region considered.

Although only FSE observations and predictions are considered in this work, the approach below applies to any other variable, related to the FSE through a linear transfer function, such as pressure, particle or surface velocity, buoy or ship motion in any degree of freedom, or any combination of the above.

Under the Gaussian assumption 1, the wave field is entirely characterised by its covariance function, defined as $\textbf{r}_{\eta \eta }(\textbf{x}_1,\textbf{x}_2,t_1,t_2) := \langle \eta (\textbf{x}_1,t_1) \eta (\textbf{x}_2,t_2) \rangle$ , where $\langle . \rangle$ denotes the mathematical expectation. Because of the stationarity and homogeneity properties of assumptions 2a and 2b, the covariance function depends only on the relative position and relative instant (lag), and therefore can be written as $ \textbf{r}_{\eta \eta }(\Delta \textbf{x},\tau )$ where $\Delta \textbf{x} = \textbf{x}_2-\textbf{x}_1$ and $\tau =t_2-t_1$ . The covariance function can be calculated from the wave spectrum through an inverse Fourier transform, the Wiener–Khinchin relation (Ochi Reference Ochi2005), which, for a directional sea, takes the following form:

(3.1) \begin{equation} \textbf{r}_{\eta \eta }(\Delta \textbf{x},\tau ) = \frac {1}{2\pi }\int \limits _{0}^{\infty }\int \limits _{\theta = 0}^{2\pi } S_{\eta \eta }(\omega ,\theta ) e^{j[\textbf{k}(\theta ,\omega ).\Delta \textbf{x} - \omega \tau ]}{\textrm{d}}\theta {\textrm{d}}\omega \end{equation}

where $\textbf{k}(\theta ,\omega )$ is a vector with direction $\theta$ and magnitude $k$ , related to $\omega$ through a dispersion relation, and $S_{\eta \eta }(\omega ,\theta )$ is the frequency-direction wave spectrum.

3.2. Optimal predictor calculation

From assumptions 1 and 2, the observed and predicted vector time series $\textbf{z}_o(t)$ and $\textbf{z}_p(t)$ , defined in § 2, jointly form a stationary, Gaussian, multi-dimensional process, characterised by the autocovariance matrix-valued functions $\textbf{r}_{oo}(\tau ):=\langle \textbf{z}_o(t)\textbf{z}_o^{\intercal} (t+\tau )\rangle$ and $\textbf{r}_{pp}(\tau ):=\langle \textbf{z}_p(t)\textbf{z}_p^{\intercal} (t+\tau )\rangle$ , and the cross-covariance function $\textbf{r}_{op}(\tau ):=\langle \textbf{z}_o(t)\textbf{z}_p^{\intercal} (t+\tau )\rangle = \textbf{r}_{po}^{\intercal} (-\tau )$ (Papoulis & Pillai Reference Papoulis and Pillai2002). Therefore, at any time $t$ , the observed values $\textbf{Z}_o(t)$ from (2.1) and the values to be predicted, $\textbf{Z}_p(t)$ from (2.2), are zero-mean, jointly Gaussian vectors, entirely characterised by their covariance matrices, denoted as $\Sigma _{oo}$ , $\Sigma _{pp}$ and $\Sigma _{op} = \Sigma _{po}^{\intercal}$ .

The covariance matrices $\Sigma _{oo}$ , $\Sigma _{pp}$ and $\Sigma _{op} = \Sigma _{po}^{\intercal}$ do not depend on $t$ , and they are block-Toeplitz, constructed as follows:

(3.2) \begin{align} &\Sigma_{oo} = \begin{pmatrix} \textbf{A}_{1,1} & \cdots & \textbf{A}_{1,N_o}\\ \vdots & \ddots & \vdots \\ \textbf{A}_{N_o,1} & \cdots & \textbf{A}_{N_o,N_o} \end{pmatrix} \in \mathbb {R}^{M_oN_o\times M_oN_o}\nonumber\\ &\Sigma _{pp} = \begin{pmatrix} \textbf{B}_{1,1} & \cdots & \textbf{B}_{1,N_p} \\ \vdots & \ddots & \vdots \\ \textbf{B}_{N_p,1} & \cdots & \textbf{B}_{N_p,N_p} \end{pmatrix} \in \mathbb {R}^{M_pN_p\times M_pN_p} \\ &\Sigma _{op} = \begin{pmatrix} \textbf{C}_{1,1} & \cdots & \textbf{C}_{1,N_p} \\ \vdots & \ddots & \vdots \\ \textbf{C}_{N_o,1} & \cdots & \textbf{C}_{N_o,N_p} \end{pmatrix} \in \mathbb {R}^{M_oN_o\times M_pN_p} \nonumber\\ &\Sigma _{po}= \Sigma _{op}^{\intercal} \nonumber\end{align}

with blocks $\textbf{A}$ , $\textbf{B}$ and $\textbf{C}$ defined as

\begin{equation} \begin{cases} \textbf{A}_{i,j} = \textbf{r}_{oo}(\tau _o(j)-\tau _o(i)), \ \forall (i,j) \in [1\ ..\ N_o]^2 \\ \textbf{B}_{i,j} = \textbf{r}_{pp}(\tau _p(j)-\tau _p(i)), \ \forall (i,j) \in [1\ ..\ N_p]^2\\ \textbf{C}_{i,j} = \textbf{r}_{op}(\tau _p(j)-\tau _o(i)), \ \forall (i,j) \in [1\ ..\ N_o]\times [1\ ..\ N_p] \end{cases}.\end{equation}

Under such assumptions, the mean value of $\textbf{Z}_p(t)$ given the available observations $\textbf{Z}_o(t)$ , denoted $\mathbf{\mu }_{p|o}(t)$ , constitutes the optimal predictor for $\textbf{Z}_p(t)$ , in a mean-square error sense, given the observations $\textbf{Z}_o(t)$ . Here $\mathbf{\mu }_{p|o}(t)$ is calculated as follows (Papoulis & Pillai Reference Papoulis and Pillai2002):

(3.3) \begin{equation} \mathbf{\mu }_{p|o}(t) = \textbf{P} \textbf{Z}_o(t) .\end{equation}

Here $\textbf{P}$ is called the optimal prediction matrix, and is calculated as follows:

(3.4) \begin{equation} \textbf{P} = \Sigma _{po}\Sigma ^{\dagger} _{oo} \end{equation}

and the $\dagger$ (dagger) superscript denotes the inverse – if $\Sigma _{oo}$ is full-rank – or a pseudoinverse – if $\Sigma _{oo}$ is not full-rank, which indicates that the information in $\textbf{Z}_o(t)$ is statistically redundant, as shown by Mérigaud & Ringwood (Reference Mérigaud and Ringwood2018a ).

The corresponding prediction error, $\textbf{E}(t):=\mathbf{\mu }_{p|o}(t)-\textbf{Z}_p(t)$ , is zero-mean Gaussian, characterised by the optimal prediction error covariance matrix $\Sigma _{p|o}$ calculated as follows:

(3.5) \begin{equation} \Sigma _{p|o} = \Sigma _{pp}-\Sigma _{po}\Sigma ^{\dagger} _{oo}\Sigma _{op} .\end{equation}

Hence, at each time step $t$ , the optimal prediction for times $(t+\underline {\tau }_p \ldots t+\overline {\tau }_p)$ is obtained through a simple multiplication, by $\textbf{P}$ , of the observed data at times $(t+\underline {\tau }_o \ldots t+\overline {\tau }_o)$ .

The only remaining challenge is the calculation of the covariance matrices (3.2) at discrete lag values corresponding to the prediction configuration, but from § 3.1, the entries of $\textbf{r}_{oo}$ , $\textbf{r}_{pp}$ and $\textbf{r}_{op}$ may be computed for any lag $\tau$ using the wave spectrum information, through (3.1). For example, entry $i,j$ of $\textbf{r}_{oo}(\tau )$ is $\textbf{r}_{oo,i,j}(\tau ) = \textbf{r}_{\eta \eta }(\textbf{x}_j-\textbf{x}_i,\tau )$ , where $\textbf{r}_{\eta \eta }$ is given in (3.1). Therefore, over short time scales below a few tens of minutes, the optimal prediction matrix $\textbf{P}$ and the optimal prediction error covariance matrix $\Sigma _{p|o}$ are both time-invariant, and only need to be updated as the wave condition (i.e. the spectrum) evolves, thus slowly modifying the statistical properties of all measured and predicted processes. The prediction error covariance matrix provides the expected prediction performance for each time horizon. Finally, although the physics of wave propagation is not explicit in the prediction (3.4), it is implicitly modelled when calculating the covariance structure, through the dispersion relation which appears in (3.1).

3.3. Accounting for current

The presence of current, or, equivalently, of a non-zero forward speed if the measurement system is mounted on a moving vessel, result in a Doppler shift, which implies that waves with a given wavelength and direction have an apparent frequency which deviates from their intrinsic frequency (i.e. their frequency in the absence of current, typically determined by the linear dispersion relation). For example, a single wave component is written as $\eta (\textbf{x},t) = \Re \{\hat {\eta } e^{j[\textbf{k}.(\textbf{x} - \textbf{U} t) - \omega t]}\}$ , where $\textbf{U}$ is the current vector, so that an observer moving at speed $\textbf{U}$ would see a ‘usual’ progressive wave $\eta (\textbf{x},t) = \Re \{\hat {\eta } e^{j[\textbf{k}.\textbf{x} - \omega t]}\}$ . Accordingly, the Wiener–Khinchine relation (3.1) must be amended as follows to account for current:

(3.6) \begin{equation} \textbf{r}_{\eta \eta }(\Delta \textbf{x},\tau ) = \frac {1}{2\pi }\int \limits _{0}^{\infty }\int \limits _{\theta = 0}^{2\pi } S_{\eta \eta }(\omega ,\theta ) e^{j[\textbf{k}(\theta ,\omega ).\Delta \textbf{x} - [\textbf{k}(\theta ,\omega ).\textbf{U} + \omega ] \tau ]}{\textrm{d}}\theta {\textrm{d}}\omega .\end{equation}

Note that in (3.6), the spectrum is the intrinsic spectrum, i.e. the spectrum which would be measured if there were no current. It can be found from (3.6) that wave measurements at different locations still form together a stationary, Gaussian multidimensional process, with covariance depending only on the lag between measurement instants and on relative positions. Only the computation of the process covariance, by using (3.6) instead of (3.1), is affected.

The above formulations (3.1) and (3.6), to obtain point-to-point covariance values, rely on the assumption of homogeneous current and wave fields. Such a model could be reasonable at scales below 10 km in the open sea Ardhuin et al. (Reference Ardhuin, Gille, Menemenlis, Rocha, Rascle, Chapron, Gula and Molemaker2017), but definitely not in intermediate and shallow depths with locally complex bathymetry or topography, which can be found in coastal areas. In the latter cases, the current, the dispersion relation, as well as the wave spectrum, could be varying within the spatial domain considered, which would considerably complicate the calculation of point-to-point covariance. In particular, the fact would no longer hold, that the spatial dependency of the covariance function $\textbf{r}_{\eta \eta }$ is only in the relative position $\Delta \textbf{x}$ , thus requiring a formulation depending on absolute positions, i.e. $\textbf{r}_{\eta \eta }(\textbf{x}_1,\textbf{x}_2,\tau )$ . Note, however, that if the covariance $\textbf{r}_{\eta \eta }(\textbf{x}_1,\textbf{x}_2,\tau )$ can still be computed, either from direct calculations, or from observations, then the rest of the statistical approach, in (3.2) to (3.6), remains unchanged.

4.1. Temporal observation interval

First, FSE measurements are taken in a single observation location $x = 0$ , for five peak wave periods into the past, i.e. $\underline {\tau }_o = -5T_p$ and $\overline {\tau }_o = 0$ , with a sampling time $T_p/10$ . The probabilistic approach of § 3 is applied for many prediction locations downwave along the $x$ axis, with a prediction time interval defined by $\underline {\tau }_p = -5T_p$ and $\overline {\tau }_p = 10T_p$ . For each prediction location, the prediction uncertainty is quantified in terms of NMSE, obtained by dividing the diagonal terms of the prediction covariance matrix $\Sigma _{p|o}$ (3.6) by $m_0$ , the variance of $\eta$ .

The NMSE map, resulting from the above procedure, is first visualised in figure 3(a) in the distance–time continuum. By drawing the NMSE contour lines for specified thresholds as in figure 3(b), the probabilistic analogues to deterministic ‘predictable zones’ (DPZs) are obtained. From this point on, PPZs shall always be defined with respect to an NMSE threshold, expressed in %, for example ‘the 10% PPZ’. Those PPZ contour lines, as can be observed in figure 3(b), form triangular zones, similarly to the deterministic theory illustrated in figure 1(a). There remains, like in DPZs, to assess whether each PPZ can be meaningfully related to two limiting group velocities, representative of the wave spectrum.

Figure 3.

Two-dimensional PPZs in the distance–time continuum, where $\tilde {x} = x/\lambda _p$ and $\tilde {\tau }=\tau /T_p$ : (a) NMSE map; (b,d) PPZs defined by NMSE contour lines. The PPZ lower slopes are retrieved from an energy criterion based on the upper frequency range: (c) let $f_{z\%}$ be the frequency above which $z\,\%$ of the spectrum energy are confined, and $c_g^{z\,\%}$ the corresponding group velocity. Then the $z\, \%$ PPZ is bounded from below by a slope $\approx 1/c_g^{z\,\%}$ , see (b) and (d). In contrast, the PPZ upper limit is poorly approximated by (b) a similar energy criterion in the lower frequency range, but well approximated by (d) the group velocity at peak frequency.

Beginning with the slowest group velocity, which determines the lower limit of each triangular predictable zone, it is found that a simple energy criterion, illustrated in figures 3(b) (dashed lines) and 3(c), provides a satisfactory agreement: considering the $z\, \%$ PPZ, the lower limiting group velocity $c_g^{z\,\%}$ is found at the frequency $f_{z\,\%}$ , above which $z\,\%$ of the spectrum integral area is confined. To the best of the authors‵ knowledge, this is the first time a quantitative criterion is proposed, to relate a particular cutoff frequency to the predictable zone.

Now examining the fastest group velocity, which should determine the upper limit of each triangular predictable zone, applying a similar energy criterion does not provide convincing results, as can be appreciated from the dashed–dotted lines in figure 3(b). In fact, the group velocities, which can be calculated from the lowest (subpeak) frequency content of the spectrum, provide exceedingly conservative estimates of the predictable zones, which is a result similar to those of Law et al. (Reference Law, Santo, Lim and Chan2020), Hlophe et al. (Reference Hlophe, Wolgamot, Kurniawan, Taylor, Orszaghova and Draper2021) and Huchet et al. (Reference Huchet, Babarit, Ducrozet, Ferrant, Gilloteaux and Droniou2023), discussed in the Introduction. In contrast, taking the peak group velocity to determine the upper limit of the triangular zones yields satisfactory agreement for all threshold values considered here, as can be appreciated in figure 3(d) (cyan, solid line). The latter result, which is novel and provides more optimistic prediction horizons, also differs slightly from those of Naaijen et al. (Reference Naaijen, Trulsen and Blondel-Couprie2014), where the limiting upper slopes at different error thresholds appear significantly different from each other.

This can be understood by going back to the line of reasoning which leads to deterministic prediction diagrams. For a wave group to be predicted, it must have travelled through the observation location $x=0$ during the observation time interval. Therefore, for a point $(x,t)$ to belong to the predictable zone, all wave groups, which together contribute to the wave signal $(x,t)$ , must have travelled through the observation location during the observed time interval. Such a train of thoughts implies that the wave process is decomposed into many wave groups, or equivalently, that the wave spectrum has been partitioned into many, sharply concentrated spectra (each of which corresponding to a random wave group). In the upper frequency region, the spectrum decreases gradually, and therefore such a partitioning can be seen as a continuum of many, small wave groups, which would explain why the criterion based on the energy integral is satisfactory. In contrast, because the spectrum decreases abruptly below the peak frequency, the spectrum peak can be seen as the spectrum of only one random wave group, centred around the peak frequency, and accounting for a large fraction of the total spectrum energy. Such an interpretation is confirmed in appendix A, by constructing a fictitious spectrum, with a more gradual decrease of the energy content below the peak frequency: In the latter case, the upper PPZ limits can also be inferred from an energy criterion in the lower frequency range, similarly to the lower PPZ limit.

4.2. Spatial observation interval

In § 4.1, the prediction is initialised using wave history at a single point. However, many RTWF methods rather rely on a spatial initialisation of the prediction, whereby a snapshot of the sea surface, at a single instant, is used as an input to the forecasting model. The very same methodology outlined in §§ 2 and 3 is followed to study such spatially distributed observations. More specifically, the observation time interval reduces to a single instant, i.e. $\underline {\tau }_o = \overline {\tau }_o = 0$ in the notations of § 2, while the list of observation locations $\textbf{x}_1, \ldots \textbf{x}_{M_o}$ spans the spatial interval of the observation snapshot. figure 4 shows the results obtained for observation points spanning two peak wavelengths, with 20 points per peak wavelength.

Figure 4.

Two-dimensional PPZs for a ‵snapshot‵ measurement configuration. Lower and upper limits of PPZs can be retrieved from the same group velocities as in the configuration of figure 3. See figure 3(c) for colour interpretation of the straight lines.

As can be appreciated in the figure, in the spatially distributed scenario, the PPZs can be approximated using the same two group velocities as in the ‘time-history’ scenario: the peak group velocity, together with the group velocity of the upper cutoff frequency, associated with the prescribed NMSE threshold.

In summary, the probabilistic approach is consistent with the two-dimensional deterministic prediction diagram, which the results of Mérigaud & Ringwood (Reference Mérigaud and Ringwood2018b ) already suggested in a qualitative way. More specifically, the present results confirm that, in a two-dimensional scenario, predictable zones can be constructed from representative group velocities in the wave spectrum, with the novel clarifications that the fastest group velocity to take into account should be the peak group velocity, and that the slowest group velocity is determined quantitatively, based on an energy criterion, depending on the threshold retained to define the predictable zone.

5.1. Single observation point

It is first examined how the waves can be predicted using a single measurement location, like in the two-dimensional case of § 4.1. Considering the spectrum with the least directional spreading ( $s_{max}=75$ ), and observations of five peak wave periods into the past, figures 5(a) and 5(b) represent the 10 % and 20 % PPZs. Also shown on the figures, are the two planes which would limit the predictable zone, if the waves were unidirectional (see figure 3 d for colour interpretation), as well as a vertical, cylindrical envelope, of which the meaning is discussed next.

Figure 5.

10 % and 20 % PPZs with $s_{max} = 75$ and a single observation point in $(x,y) = (0,0)$ . Waves are observed for five peak periods into the past. Vertical spacing between $(\tilde {x},\tilde {y})$ contour lines: $\Delta \tilde {\tau } = 1/5$ . The vertical cylindrical domains, generated using 10 % and 20 % contour lines of figure 6(c), are also represented. (a)10% probabilistic predictable zones. Upper and lower limiting slopes of the two-dimensional case are represented by planes with slopes $1/c_g^{p}$ and $1/c_g^{10\% \text{upper}}$ , respectively. See Figure 3c for color interpretation (b) 20\% probabilistic predictable zones. Upper and lower limiting slopes of the two-dimensional case are represented by planes with slopes $1/c_g^{p}$ and $1/c_g^{20\% \text{upper}}$ , respectively. See Figure 3c for color interpretation.

As can be appreciated in figure 5, the spatial and temporal extent of the predictable zone is dramatically reduced, with respect to the unidirectional case. Admittedly, at very short distances from the observation point (less than 0.5 $\lambda _p$ in this case), the upper and lower limits of the predictable zone are well approximated by the two limiting planes. However, the predictable zone rapidly shrinks with the distance from the prediction point, so that the maximum time horizon for the 10 % NMSE criterion is smaller than $T_p$ . By loosening the prediction quality criterion to 20 % (figure 5 b), the predictable zone is expanded, but still collapses to zero more rapidly than the two limiting planes would predict.

It is thus clear from figure 5 that directional spreading has a strong adverse effect on wave predictability, when a single measurement location is employed. In order to better isolate the effect of directional spreading from that of wave dispersion, we can consider the following hypothetical scenario, whereby waves are predicted at a given instant (without loss of generality, the present instant $\tau =0$ ), at a given location $(x_p,y_p)$ of the horizontal plane, using waves observed in $x=y=0$ measured both in the past and in the future, as sketched in figure 6(a). Furthermore, we can make the observation duration grow to infinity, i.e. $\underline {\tau }_o \rightarrow -\infty$ and $\overline {\tau }_o \rightarrow \infty$ . In practice, $\overline {\tau }_o$ and $\underline {\tau }_o$ are set to large enough, but finite values $\overline {\tau }_o = -\underline {\tau }_o = 20$ $T_p$ : beyond $\overline {\tau }_o = -\underline {\tau }_o = 10$ $T_p$ there is little or no change in the results, within the spatial domain considered. In other words, we quantify how well the wave signal at a location $(x_p,y_p)$ can be reconstructed, knowing the whole time history, both past and future, of the wave signal in $(0,0)$ . We also hereby suppress the time dependency of the problem.

By making $(x_p,y_p)$ vary across the horizontal plane and calculating the corresponding NMSE, the contour lines of figures 6(b) and 6(c) are obtained. Note that, because the results are symmetric with respect to the $y$ axis, only positive $x$ values are shown.

Figure 6.

Schematic diagram of wave prediction at the present time in $(x_p,y_p)$ , using both past and future observations in $x=y=0$ (figure 6 a); corresponding achievable NMSE contour lines in the horizontal plane when the observation time interval grows to infinity, for two different spreading parameters $s_{max}$ (figures 6 b and 6 c). (b) NMSE contour lines, $s_{max} = 10$ (c) NMSE contour lines, $s_{max} = 75$ .

figure 6 shows that the achievable prediction performance decreases with the distance to the observation location, most rapidly in the direction transverse to the main wave direction, and least rapidly in the main wave direction. The decrease is also stronger for widely spread seas (figure 6 b) than for long-crested seas (figure 6 c). Clearly, the prediction NMSE in figure 6, obtained with an infinitely long observation time interval, should be smaller than the NMSE obtained with observations, limited to a finite time interval. Therefore, each contour line in figures 6(b) and 6(c) defines a spatial domain, which limits the extent of the corresponding PPZ for a single observation location, regardless of the observed time segment.

To visualise the latter statement, figures 5(a) and 5(b) also show a representation of the 10 % (respectively, 20 %) contour lines of figure 6(c), extended vertically to form a cylindrical domain. As evidenced in the figure, the three-dimensional predictable zone for a single observation location can be constructed from three ingredients: two limiting group velocities, as in the unidirectional case, and a cylindrical envelope, centred around the observation location, which limits the PPZ spatial extent under the effect of directional spreading.

Considering that the results of figure 5 are obtained in the most favourable scenario $s_{max} = 75$ , it seems unlikely that forecasts of any practical interest can be achieved beyond one wave period into the future, when only one observation location is used. Another striking feature of those results is that they do not resemble the fan-shaped, DPZs, exemplified in figure 1(b). It is therefore important to understand the difference between the results of the present theory, and those of the deterministic approach.

5.2. Comparison with the deterministic theory

For the single-observation location scenario, DPZs calculated following Qi et al. (Reference Qi, Wu, Liu and Yue2018b ) are compared with the PPZs of the present modelling framework. A first challenge is that, in the present work, unlike in the deterministic approach, PPZs are defined with respect to a prescribed NMSE level. In the following examples, without loss of generality, a 10 % NMSE threshold is chosen as a reasonable approximation to a ‘perfect’ prediction. Another challenge is to choose the two limiting frequencies of the deterministic approach, so as to make the deterministic and probabilistic cases as closely comparable as possible. To that end, the deterministic approach is implemented using the peak frequency $f_p$ as the lower cutoff, while the upper cutoff frequency is set to $f_{10\%}$ , accordingly to the results of § 4. This ensures that the deterministic and PPZs will be limited by the same two group velocities. Finally, the three-dimensional deterministic theory requires that all the spectrum energy propagates between two cutoff directions, $\pm \Delta _{\theta }$ , where $\Delta _{\theta } \leqslant \pi /2$ . Thus, for a fair comparison, the probabilistic results are also calculated based on a directionally truncated spectrum $S_{\text {trunc}}(f,\theta )$ . More specifically, the JONSWAP–Mitsuyasu directional spectrum $S(f,\theta )$ of § 5.1 (with the most optimistic parameter $s_{max} = 75$ ) is directionally truncated,

(5.1) \begin{equation} S_{\text {trunc}}(f,\theta ) = \begin{cases} S(f,\theta ) & |\theta | \lt \Delta _{\theta } \\ 0 & |\theta | \geqslant \Delta _{\theta } \end{cases} .\end{equation}

The value $\Delta _{\theta } = \pi /6$ is chosen in the calculation of both PPZs and DPZs, and two different observation duration values 5 and 10 $T_p$ are considered. The resulting 10 % PPZs are shown in figures 7(a) and 7(c), while the corresponding DPZs are shown in figures 7(b) and 7(d). The two limiting planes with slopes $1/c^p_g$ and $1/c^{10\,\% \text {upper}}_g$ , which bound the predictable zones in both cases, are also represented. It is immediately clear that the DPZs are not only different in shape, but also significantly larger in size than their analogues defined in a probabilistic sense. Using five peak wave periods of observations, the deterministic theory predicts that waves can be forecast accurately up to $3.5$ peak wave periods into the future (figure 7 b), versus only $1.5$ following the probabilistic theory (figure 7 a). Using 10 peak wave periods of observations, instead of five, doubles the achievable forecasting horizon in the deterministic theory (figure 7 d), while it leaves it unchanged in the probabilistic theory (figure 7 c).

Figure 7.

10 % PPZs with $s_{max} = 75$ and a single observation point in $(x,y) = (0,0)$ . Vertical spacing between the contour lines: $\Delta \tilde {\tau } = 0.5$ (0.5 $T_p$ ). The directional wave spectrum is truncated for $|\theta |\geqslant \pi /6$ . Upper and lower limiting slopes of the two-dimensional case are represented by planes with slopes $1/c_g^{p}$ and $1/c_g^{10\,\% \text {upper}}$ , respectively. See figure 3(c) for colour interpretation. (a) 10% probabilistic predictable zones. Waves are observed for 5 peak periods into the past (b) Deterministic predictable zones. Waves are observed for 5 peak periods into the past (c) 10% probabilistic predictable zones. Waves are observed for 10 peak periods into the past (d) Deterministic predictable zones. Waves are observed for 10 peak periods into the past.

In fact, it is an inherent property of the DPZ construction method, that the achievable time horizon $T_{max}$ is proportional to the observation duration $T$ . The relation is

(5.2) \begin{equation} T_{max} = \frac {T}{\frac {c^{{max}}_g}{c^{{min}}_g \cos {\Delta _{\theta }}}-1} \end{equation}

which reduces to $T_{max} = {T}/({c^{{max}}_g/c^{{min}}_g-1})$ in the unidirectional case (Qi et al. Reference Qi, Wu, Liu and Yue2018b ; Kim et al. Reference Kim, Ducrozet, Bonnefoy, Leroy and Perignon2023). Therefore, by increasing the observation duration at a single location, the forecasting horizon could, in principle, be increased to arbitrarily large values, at a downwave location farther and farther away from the observation point. The latter property is in clear contradiction with the probabilistic theory, in which predictable zones are spatially bounded under the effect of directional spreading, following the error maps shown in figure 6.

In the deterministic theory, the construction of the three-dimensional predictable zone is based on the intersection of two-dimensional predictable zones, with incident wave directions spanning the interval $[-\Delta _{\theta }, \Delta _{\theta }]$ . However, for this construction to hold in a real-sea scenario, the forecasting model should be able to decompose, within the observed wave signal, the contributions from each wave direction, before propagating them in time and space. In a realistic sea, those contributions are randomly mixed at the observation location, so that it is not possible to perfectly decompose the wave signal into its directional components. Such a ‘directional mixing’ effect is, implicitly, modelled in the probabilistic theory, but not in the deterministic representation.

Therefore, in the present probabilistic model, with a single observation location, increasing the forecasting horizon cannot be achieved by looking further into the past. Instead, resolving for the different incoming wave directions can only be achieved by spatially extending the domain of observation, which is the topic of the next section.

5.3. A line of observation points

Probabilistic predictable zones are now examined for a straight line of observation points, extending in the direction transverse to the peak wave direction $\theta _p$ . The line length is set to two peak wavelengths. Besides, it is found that no appreciable change can be observed in prediction error beyond a point density of 10 points per peak wavelength, which is therefore chosen throughout this section. The waves are observed for five peak wave periods into the past, with a time step of $T_p$ /10.

figure 8 shows the 10 % and 20 % PPZs, for two different values of the spreading parameter $s_{max}$ , along with the two limiting planes corresponding to the two-dimensional case. In comparison with the single-point configuration, the spatial information provided by the line of measurements overcomes the adverse effect of directional mixing, allowing prediction horizons ranging from 1.5 $T_p$ (in the least favourable case $s_{max} = 10$ ) to 4 $T_p$ ( $s_{max} = 75$ ), for a 10 % NMSE criterion.

Figure 8.

10 % and 20 % PPZs with a line of observation points. Waves are measured for five peak wave periods into the past. Vertical spacing between the contour lines: $\Delta \tau = 0.5 T_p$ . Upper and lower limiting slopes of the two-dimensional case are represented by planes with slopes $1/c_g^{p}$ and $1/c_g^{10\,\% \text {upper}}$ , respectively. The effect of directional spreading is approximated by cylindrical envelopes (vertically stretched 10% and 20% GOP contour lines). See figures 3(c) and 9 for colour interpretation. (a) 10 % PPZ, $s_{max} = 10$ (b) 20 % PPZ, $s_{max} = 10$ (c) 10 % PPZ, $s_{max} = 75$ (d) 20 % PPZ, $s_{max} = 75$ .

Yet, along with the two limiting wave groups, directional spreading maintains a shrinking effect onto the predictable zone geometry, although at a significantly slower rate than in the single-point case. In the same way that wave group travelling speeds govern upper and lower limits of predictable zones, it is interesting to outline a similar physical interpretation for the directional shrinking effect, based on wave group travelling directions.

To that end, we shall make the simplifying assumption that the wave field is made of random wave groups, centred around the spectrum peak frequency, and idealised as points representing the wave group centre. The group travelling directions are randomly drawn in $[0; 2\pi ]$ , with a probability distribution represented by the directional spreading function. We also approximate the discrete measurement points as a continuous line, stretching partway in the transverse direction across a length $L_y$ , as sketched in figure 10.

Figure 9.

Geometrical calculation of the group observation probability for a point $P(x,y)$ .

Next, we consider one such wave group, arriving at a point $P$ of coordinates $(x,y)$ at $t=0$ . A necessary condition, for this wave group to be predicted in $(x,y)$ , is that it has travelled through the line of observations some time in the past. The probability that it is the case can be calculated geometrically, from the directional distribution function at the peak frequency, as pictured in figure 10. With the notations of figure 10, the ‘group observation probability’ (GOP) is mathematically expressed as follows:

(5.3) \begin{equation} \text {GOP}(x,y) = \int \limits _{\mathcal {S}(x,y)} D_p(\theta ) {\textrm{d}}\theta \end{equation}

where $D_p(\theta )$ is the directional spreading function at the peak frequency, and $\mathcal {S}(x,y)$ is the angular sector comprised between the two extreme wave directions traversing both the observation line and the point $P$ of coordinates $(x,y)$ . Repeating the calculation for all points $(x,y)$ in the horizontal plane, we obtain a GOP map, representing the fraction of signal energy which has travelled through the observation line prior to reaching each point. The resulting contour lines, for selected GOP thresholds, can be visualised in figures 9(a) and 9(b) using normalised coordinates $\tilde {x} = x/\lambda _p$ and $\tilde {y} = y/\lambda _p$ .

Figure 10.

Group observation probability $\text {GOP}(\tilde {x},\tilde {y})$ , i.e. fraction of directional energy crossing the observation line before reaching point $(\tilde {x},\tilde {y})$ , for a strip of measurements extending indefinitely in the time direction. The directional spreading function is taken at the peak frequency. (a) GOP contour lines, $s_{max} = 10$ (b) GOP contour lines, $s_{max} = 75$ .

Stretching those GOP contour lines vertically, as seen in figure 8, provides, in combination with the two limiting planes, a reasonable approximation of the predictable zone satisfying the corresponding thresholds. Thus, although the GOP map certainly relies on strong approximations, it is straightforward to calculate, and it seems sufficient to provide an estimate of the spatial bounds which limit PPZs under the effect of directional spreading.

5.4. A square snapshot of observation points

As mentioned in § 4.2, many wave prediction methods rely on a snapshot of the sea surface to initialise the prediction, rather than on a long time history. Observations taken in a square area of the sea surface are considered here. The square side length is set to one peak wavelength. As in the previous section, the point density in the transverse direction is set to 10 points per peak wavelength, while in the longitudinal direction, a higher density of 20 points per peak wavelength is found preferable to achieve convergence of results.

Figure 11 shows the 10 % PPZs, for two different values of the spreading parameter $s_{max}$ , along with the two limiting planes corresponding to the two-dimensional case. While being transported, and gradually shrinking in the longitudinal direction $x$ between the two limiting planes, the PPZ also shrinks in the transverse direction $y$ , thus taking a pyramidal aspect. The transverse shrinking rate is faster for short-crested waves, and more gradual for long-crested waves, as can be appreciated by comparing figures 11(a) and 11(b).

Figure 11.

10 % PPZs with a line of observation points. Waves are measured for five peak wave periods into the past. Vertical spacing between the contour lines: $\Delta \tau = 0.5 T_p$ . Upper and lower limiting slopes of the two-dimensional case are represented by planes with slopes $1/c_g^{p}$ and $1/c_g^{10\,\% \text {upper}}$ , respectively. See figures 3(c) and 12 for colour interpretation. (a) 10 % PPZ, $s_{max} = 10$ (b) 10 % PPZ, $s_{max} = 75$ .

Similarly to the line arrangement, it is desirable to determine a simple and meaningful approximation to the directional spreading effect. To that end, we shall make the same simplifying assumptions of random wave groups, centred around the spectrum peak frequency, and idealised as points travelling at the peak group velocity. The travelling direction is randomly drawn in $[0; 2\pi ]$ , with a probability distribution represented by the directional spreading function. In addition, the domain of observation is simplified to a strip, extending indefinitely in the longitudinal direction, hence suppressing the $x$ dependency of the problem, as sketched in figure 13.

Figure 12.

Geometrical calculation of the GOP for a point $P(y,\tau )$ .

Figure 13.

Group observation probability $\text {GOP}(\tilde {y},\tilde {\tau })$ , i.e. fraction of directional energy crossing the strip of measurement at time $\tilde {\tau } = 0$ before reaching point $(\tilde {y},\tilde {\tau })$ , for a strip of measurements extending indefinitely in the $x$ direction. The directional spreading function is taken at the peak frequency. (a) GOP contour lines, $s_{max} = 10$ (b) GOP contour lines, $s_{max} = 75$ .

We then examine a wave group reaching a point $P$ of coordinates $(y,\tau )$ , and we calculate its observation probability, i.e. the probability that it has crossed the observation domain at the measurement instant $\tau =0$ . Since the wave group travels at speed $c_g^p$ , we know that, at $\tau =0$ , the group was located along the circle of radius $\tau c_g^p$ , centred in $y$ . The wave group has thus been observed, if and only if its travelling direction belongs to the angular sector $\mathcal {S}$ , of the circle intersection with the observation strip. The probability that it is so is calculated by integrating the directional spreading function, restricted to $\mathcal {S}$ . Using the notations of figure 13, this is mathematically expressed as follows:

(5.4) \begin{equation} \text {GOP}(y,\tau ) = \int \limits _{\mathcal {S}(y,\tau )} D_p(\theta ) {\textrm{d}}\theta \end{equation}

with

(5.5) \begin{equation} \mathcal {S}(y,\tau ) = \{\theta : |\Re \{y + \tau c_g^p e^{j\theta }\}| \leqslant \frac {L_y}{2} \} .\end{equation}

Hence, similarly to the line configuration, a two-dimensional map $\text {GOP}(y,\tau )$ is obtained, representing, in each point $(y,\tau )$ , the fraction of signal energy which has been observed at the measurement instant. Drawing the corresponding contour lines for selected thresholds (figures 12 a and 12 b), and stretching them horizontally in the $x$ direction, cylindrical domains aligned with the $x$ axis are obtained, which are visualised in figure 11. In combination with the two limiting planes, they describe the gradual shrinking of the predictable zones with reasonable accuracy.

Like in the line observation configuration, the PPZ geometry can thus be interpreted and approximated through relatively simple physical considerations. Lines and squares are obviously not the only possible types of measurement configurations, but they are archetypal examples, which can certainly provide guidelines to determine the necessary measurement ranges in time and in the two horizontal dimensions, depending on the prediction horizon and accuracy requirements.

6.1. Dataset presentation

The dataset employed in this study consists of a three-dimensional movie of a wave field, reconstructed by means of a system of two stereo cameras, from an ocean observation platform in the Yellow Sea in Korea, and thoroughly described and analysed by Benetazzo et al. (Reference Benetazzo, Bergamasco, Yoo, Cavaleri, Kim, Bertotti, Barbariol and Shim2018). It is part of a larger ensemble of datasets, introduced by Veras Guimãraes et al. (Reference Veras Guimãraes, Ardhuin, Bergamasco, Leckler, Filipot, Shim, Dulov and Benetazzo2020), recorded at various locations worldwide and made freely available online. Unlike other datasets in (Veras Guimãraes et al. Reference Veras Guimãraes, Ardhuin, Bergamasco, Leckler, Filipot, Shim, Dulov and Benetazzo2020), the Yellow Sea record has a spatial coverage, in the order of two peak wavelengths, large enough for relatively accurate wave spectrum estimation through direct Fourier analysis. Furthermore, both the spatial coverage and time duration allow the empirical assessment of predictable zones. The statistical analysis reported by Benetazzo et al. (Reference Benetazzo, Bergamasco, Yoo, Cavaleri, Kim, Bertotti, Barbariol and Shim2018) shows skewness and kurtosis values of $m_3=0.14$ and $m_4=3.10$ , respectively, indicative of detectable deviation from Gaussianity due to bound harmonics. However, as shall be observed in the rest of this section, such a deviation does not appear to have an appreciable negative impact on the applicability of the – albeit linear – probabilistic approach of § 3. Overall, this dataset is an almost-ideal playing field to assess wave predictability in a real sea scenario.

The sea state is characterised by a significant wave height $H_{m_0}$ of approximately 2 m, and a peak wave period $T_p$ below 6 s. In addition, as documented by Benetazzo et al. (Reference Benetazzo, Bergamasco, Yoo, Cavaleri, Kim, Bertotti, Barbariol and Shim2018), a strong current is present, with speed larger than 1 m s^–1. The water depth is approximately 50 m, which is close to deep-water conditions for the wavelengths present in the wave spectrum. Hence, bathymetry variations should have little or no impact on experimental results.

The original record is 440 s long, with a time step of 0.1 s. For the remainder of this study, a ${120}\,\textrm {m} \times {120}\,\textrm {m}$ square area is extracted, and resampled regularly in space with a 2 m resolution in both $x$ and $y$ directions. The wave elevation is also resampled with an 0.25 s time step. figure 14 shows a sample snapshot of the record, in the square subdomain extracted for the present study.

Figure 14.

A $120\text {m} \times 120\text {m}$ snapshot of the FSE record extracted from the dataset of Veras Guimãraes et al. (Reference Veras Guimãraes, Ardhuin, Bergamasco, Leckler, Filipot, Shim, Dulov and Benetazzo2020). The intrinsic peak wave direction and current direction are also represented.

The record is first used to estimate the underlying wave spectrum, which is accomplished using Fourier transform in the three dimensions ( $x$ , $y$ , $t$ ). One challenge is that, for the calculations of (3.6), the intrinsic spectrum must be used, while the direct, three-dimensional Fourier transform of the dataset provides the apparent spectrum (i.e. modified by surface current). Therefore, the current must be estimated, which allows the correction of the apparent spectrum, as detailed in appendix C. The resulting intrinsic spectrum is shown in figure 15.

Figure 15.

Intrinsic directional spectrum estimated from the Yellow Sea dataset. (a) Directional spectrum (b) Frequency spectrum (c) Directional spreading function (DSF), frequency-averaged and at peak frequency.

The peak wave period is found to be $T_p = {5.75}\,\textrm {s}$ and the peak wave direction is $220^{\circ}$ (in the reference frame of figure 14). With a mean water depth of 50 m, the peak wavelength is 51.6 m, so that the square area covers a little more than two peak wavelengths. Finally, wave energy is widely spread around the peak wave direction. The current magnitude is 1.09 m s^–1, with a direction of $302^{\circ}$ .

6.2. Experimental comparison procedure

Using the notations of § 2, consider the observation locations $\textbf{x}_{i = 1\ldots M_o}$ and time interval $[\underline {\tau }_o, \overline {\tau }_o]$ . Let the waves be predicted at a single location $\textbf{x}^{\prime} = (x^{\prime},y^{\prime})$ , over the time interval $[\underline {\tau }_p, \overline {\tau }_p]$ . Based on the intrinsic spectrum and current, and using the Wiener–Khinchine relation (3.6), the procedure of § 3 is used to obtain the theoretical NMSE from (3.5), exactly like in §§ 4 and 5, but also the prediction matrix $\textbf{P}$ (3.4).

The time series $\textbf{z}_o(t) = (\eta (x_1,y_1,t), \ldots , \eta (x_{M_o},x_{M_o},t))^T$ and $\textbf{z}_p(t) = \eta (x^{\prime},y^{\prime},t)$ are extracted from the empirical record at the observation and prediction locations $\textbf{x}_{i = 1\ldots M_o}$ and $\textbf{x}^{\prime}$ , and the prediction is carried out in the sliding-window fashion sketched in figure 2, using (3.4) with the prediction matrix $\textbf{P}$ . The predicted value of $\eta$ at time horizon $\tau$ , calculated from the observations available at time $t$ , is denoted as $\eta _{p|o(t)}(x^{\prime},y^{\prime},t + \tau )$ . The empirical NMSE $\tilde {\epsilon }_{\text {emp}}$ at the prediction location $\textbf{x}^{\prime} = (x^{\prime},y^{\prime})$ is calculated, for each prediction horizon $\tau$ in $[\underline {\tau }_p, \overline {\tau }_p]$ , as follows:

(6.1) \begin{equation} \tilde {\epsilon }_{{emp}}(x^{\prime},y^{\prime},\tau ) = \frac {1}{m_0 [i_{{max}}-i_{{min}}+1]} \sum \limits _{i = i_{{min}}}^{i_{{max}}} [\eta _{p|o(t_i)}(x^{\prime},y^{\prime},t_i + \tau )-\eta (x^{\prime},y^{\prime},t_i + \tau )]^2, \end{equation}

where

(6.2) \begin{equation} m_0 = \frac {1}{N_t N_x N_y} \sum \limits _{i=1}^{N_t} \sum \limits _{j=1}^{N_x} \sum \limits _{k=1}^{N_y} \eta (x_j,y_k,t_i)^2 \end{equation}

is the variance of $\eta$ , and $i_{min}$ and $i_{max}$ are the first and last time indices where the prediction can be carried out, depending on the values of $\underline {\tau }_o$ , $\overline {\tau }_o$ , $\underline {\tau }_p$ and $\overline {\tau }_p$ . Allowing the prediction location $(x^{\prime},y^{\prime})$ to vary across the whole experimental domain, empirical NMSE maps are constructed, which can be directly compared with the theoretical NMSE.

In the procedure outlined above, a significant implementation detail must be considered. Indeed, the actual wave time series may be subject to measurement noise, or unmodelled nonlinear dynamics. Furthermore, the intrinsic spectrum and the current may not have been estimated perfectly. Those effects all imply that actual wave statistics can depart slightly from those predicted in the Gaussian approximation. As a consequence, some robustness must be introduced in the prediction matrix calculation, to ensure that empirical results are not too sensitive to even minor deviations from the expected statistics. To that aim, it is assumed that some level of uncorrelated, Gaussian noise is also present in the observations, which amounts to adding a positive, constant term $\sigma ^2_{\text {n}}$ to the diagonal elements of $\Sigma _{oo}$ , where $\sigma ^2_{\text {n}}$ represents the white noise variance. Without noise modelling ( $\sigma ^2_n = 0$ ), the empirical NMSE is far below the theoretical one. By gradually increasing $\sigma ^2_n$ , the empirical NMSE increases significantly, while the theoretical NMSE remains almost unchanged. Here $\sigma ^2_n$ can thus be tuned so that theoretical and empirical results closely match, as done by Mérigaud et al. (Reference Mérigaud, Herterich, Flanagan, Ringwood and Dias2018). A value of $\sigma ^2_{\text {n}} = 0.025 m_0$ (in other words, 2.5 % of the wave signal variance) is adopted throughout this experimental study, and found sufficient to ensure robustness of results.

6.3. Single observation point

Starting with a single measurement point (chosen in the centre of the observation frame), and setting $\underline {\tau }_o = -{30}\,\textrm {s}$ , $\overline {\tau }_o = {0}\,\textrm {s}$ , with a time step of ${0.5}\,\textrm {s}$ , theoretical and empirical NMSE maps in the $(x,y)$ plane are constructed. The resulting theoretical and empirical PPZs are compared in figure 16, for a 20 % threshold value, showing a very good agreement between the theoretically predicted PPZ and the empirical one. Like in § 5.1, predictable zones at specific time horizons resemble ellipses, with their principal axis almost aligned with the peak wave direction. However, because of surface current which is almost perpendicular to the peak wave direction, the PPZ orientation is slightly deviated. Regardless of the precise shape of PPZs, empirical results tend to confirm that, using a single observation location in a directional wave field, no forecasts of practical interest can be obtained, unless the prediction requirements are extremely low, in terms of time horizon and prediction quality. Most importantly, neither the theoretical, nor the empirical results, support the existence of fan-shaped predictable zones, as constructed in the deterministic theory.

Figure 16.

Theoretical and empirical PPZs for a 20 % NMSE threshold. Waves are observed at a single location, looking for 30 s into the past. The peak wave direction and current direction are represented by purple and green arrows, respectively. (a) Theoretical PPZ computed from the wave spectrum. (b) Empirical PPZ obtained by applying the forecasting model to the experimental data. Vertical spacing between $(x,y)$ contour lines: 1 s.

Figure 17.

Theoretical and empirical PPZs for 10 % and 20 % NMSE thresholds. Waves are observed along a line perpendicular to the main intrinsic wave direction, looking for 30 s into the past. (a,c) Theoretical PPZ computed from the wave spectrum. (b,d) Empirical PPZ obtained by applying the forecasting model to the experimental data. Vertical spacing between $(x,y)$ contour lines: 1 s. The effect of directional spreading is approximated by cylindrical envelopes (vertically stretched 10 % and 20 % GOP contour lines). (a,b) 10% PPZ. (c,d) 20% PPZ.

figure 16 is only one example, since only one NMSE threshold is considered, but a more complete set of results, including a comparison of theoretical and empirical NMSE maps at selected prediction horizons, is available as Supplementary materials. Those additional results all confirm the agreement between theoretical and empirical prediction performance metrics.

6.4. A line of observation points

A line of measurement points is also investigated, oriented perpendicular to the peak intrinsic wave direction, similarly to § 5.3. The line length is two peak wavelengths like in § 5.3, with a point density of 20 points per wavelength. Like in the single-point case, observations are made between $\underline {\tau }_o = -{30}\,\textrm {s}$ and $\overline {\tau }_o = {0}\,\textrm {s}$ , with a time step of 0.5 s. Theoretical and empirical results are shown in figure 17, in the form of the 10 % and 20 % PPZs. Considering the figure, the line of observations does not allow significant prediction horizons if a 10 % NMSE criterion is retained, which can be attributed to the large directional spreading of the wave energy content. Loosening the criterion to 20 % (figure 17 b), the prediction horizon is extended to almost 10 s. Like in § 5.3, cylindrical domains, obtained from GOP contour lines, are also represented, as an approximation to the effect of directional spreading. In this case, the effect of current is taken into account to calculate the GOP, as documented in appendix B. As can be appreciated in figure 17(a), the GOP 10 % contour line is not a satisfactory approximation to the spatial extent of the 10 % PPZs, both theoretical and empirical. The agreement is significantly more convincing for the 20 % PPZs in figure 17(b). Such discrepancies may be due to the fact that GOP is calculated from the directional energy content at the peak frequency, while the actual directional energy content varies significantly across frequencies, as can be appreciated in figure 15.

Figure 18.

NMSE contour lines defining PPZs for selected threshold values, along with limiting slopes determined in § 4. Waves are observed like in figure 17. Results are taken along a segment aligned with the intrinsic peak wave direction, and crossing the middle observation point. $x^{\prime}$ denotes the longitudinal distance along the segment (i.e. the distance from the observation line).

In addition to directional spreading, it is interesting to verify the two other ingredients governing wave predictability, highlighted in § 5, namely the two limiting group velocities. To that end, the NMSE contour lines are visualised along a segment, aligned with the peak wave direction and crossing the middle observation point. This results in NMSE contour lines in the $x^{\prime}-\tau$ plane (because the current is oriented almost perpendicularly to the waves, it is not expected to have a significant effect in this cross-section of the results), where $x^{\prime}$ is the distance from the middle observation point. Those contour lines are represented in figure 18 for specific NMSE threshold values. In the same figure, straight lines are also represented, with slopes corresponding to the limiting group velocities identified in § 4. At small distances, each $z\,\%$ PPZ is bounded from above by the group velocity at peak frequency, and from below by $c_g^{z\,\%, \text {upper}}$ , determined from the energy criterion of § 4. However, instead of remaining guided between the two limiting slopes, the PPZs are spatially limited under the effect of directional spreading, long before the two slopes intersect. Compared with the two-dimensional results of § 4, it is more difficult, from figure 18 alone, to assess exactly the accuracy of the proposed limiting slopes $c_g^{z\,\%, \text {upper}}$ and $c_g^{p}$ , due to the fact that the PPZs collapse rapidly under the directional spreading effect. Nevertheless, those two group velocities remain, from this set of results, very plausible approximations.

Overall, there is a good agreement between theoretical and empirical results, which is further supported by additional figures, available as Supplementary materials, including theoretical and empirical NMSE maps in the $(x,y)$ plane for selected prediction horizons.

6.5. A square snapshot of observation points

The square snapshot configuration of § 5.4 is now considered, where the square is aligned with the peak intrinsic wave direction, the square side length is one peak wavelength, and the point density in both longitudinal and transverse directions is 20 points per peak wavelength. The NMSE results are presented as a PPZ for the 20 % NMSE criterion, in figure 19. Overall, empirical results ( figure 19 b) follow those, calculated from the wave spectrum (figure 19 a). The empirical prediction performance is, however, slightly better than expected, without calling into question the overall shape and behaviour of the PPZ. Like in the line configuration, the area of better prediction accuracy is transported along the main wave direction at the peak group velocity, and in the current direction at the current speed (see appendix B for more detail). Overall, retaining a 20 % NMSE criterion, this square configuration allows forecasts up to approximately 10 s ( $\approx 1.5 T_p$ ) into the future.

Figure 19.

Theoretical and empirical PPZs for a 20 % NMSE threshold. Waves are observed in $\tau = 0$ in a $21 \times 21$ -point square, aligned with the main intrinsic wave direction, with the sides equal to one peak wavelength. (a) Theoretical PPZ computed from the wave spectrum. (b) Empirical PPZ obtained by applying the forecasting model to the experimental data. Vertical spacing between $(x,y)$ contour lines: 1 s.

Finally, figure 20 shows an example of the time evolution of the wave field, predicted in the square snapshot configuration. For selected prediction horizons, figure20(a,c) represent the actual FSE, in and around the initial observation domain, while figure 20(b,d) shows the FSE, predicted in the same area, using the square snapshot (bounded by dotted lines) taken at $\tau = 0$ . The theoretical 20 % PPZ (i.e. figure 19 a) is also materialised by means of dashed, green lines, for each sample prediction horizon. Note that, at time $\tau =0$ , the 20 % PPZ includes both the square observed snapshot, as well as the surrounding area where the FSE can be inferred from the snapshot with NMSE smaller than 20 %.

Figure 20.

An example of wave field prediction from a square snapshot of observations, at different prediction horizons $\tau$ . The square snapshot is represented by dotted lines at $\tau = 0$ . (a,c) True wave field, (b,d) predicted wave field. Green, dashed lines: theoretical 20 % PPZ.

figure 20 illustrates how the achievable prediction quality degrades with the time horizon, increasingly resembling a low pass-filtered version of the true wave field, whereby only larger wavelengths are preserved – an observation already suggested by Mérigaud & Tona (Reference Mérigaud and Tona2022) by analysing the frequency content of the prediction error signal. In addition, incoming waves mostly propagate from the top-right corner of the domain. As a result, for larger time horizons (beyond $T_p$ ), the top-right part of the wave field essentially consists of newly entered, independent wave groups, which could by no means be inferred from the observation square at initial time $\tau = 0$ . Since there is virtually no information about the FSE in that area of the domain for larger time horizons, the optimal forecast approaches zero, i.e. the unconditional mean of the FSE field.