# Development of Phase-Resolved Real-Time Wave Forecasting With Unidirectional and Multidirectional Seas

In-Chul Kim, Guillaume Ducrozet, Yves Perignon

► To cite this version:

In-Chul Kim, Guillaume Ducrozet, Yves Perignon. Development of Phase-Resolved Real-Time Wave Forecasting With Unidirectional and Multidirectional Seas. ASME 2023 42nd International Conference on Ocean, Offshore and Arctic Engineering, Jun 2023, Melbourne, France. 10.1115/OMAE2023-104857. hal-04490911

HAL Id: hal-04490911

https://hal.science/hal-04490911v1

Submitted on 5 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Proceedings of the ASME 2023 42nd International
Conference on Ocean, Offshore and Arctic Engineering
OMAE2023

June 11-16, 2023, Melbourne, Australia

OMAE2023-104857

DEVELOPMENT OF PHASE-RESOLVED REAL-TIME WAVE FORECASTING WITH UNIDIRECTIONAL AND MULTIDIRECTIONAL SEAS

I.-C. Kim

Nantes Université, École Centrale
Nantes, CNRS, LHEEA, UMR 6598,
F-44000 Nantes, France

G. Ducrozet

Nantes Université, École Centrale
Nantes, CNRS, LHEEA, UMR 6598,
F-44000 Nantes, France

Y. Perignon

Nantes Université, École Centrale
Nantes, CNRS, LHEEA, UMR 6598,
F-44000 Nantes, France

ABSTRACT

The development of phase-resolved real-time wave forecast-
ing is outlined. This framework is an enhancement over previous
work in that the algorithm of real-time wave prediction is extended
into multidirectional seas by including the wave measurements
and components in direction. However, the computations with
multidirectional seas become much more numerically expensive,
and hence it is often not possible to accomplish a real-time sys-
tem of nonlinear ocean wave prediction. Accordingly, we sug-
gest an improved assimilation procedure in the process of wave
reconstruction, which is proven to alleviate the computational
costs and establish the numerical stability of the Lagrangian ap-
proach. In addition, given an observation zone recorded by an
optical sensor mounted on a fixed offshore structure, we provide a
spatio-temporal prediction zone where it is suitable to obtain the
prediction of the wave field by evolving the reconstructed wave
information in time and space. In order to validate the phase-
resolved wave forecasting, we conducted a tank-scale experi-
mental campaign with unidirectional seas (long-crested waves)
and multidirectional seas (short-crested waves). Through the
comparison of model performance against the laboratory data
between unidirectional and multidirectional seas, it is confirmed
that the directional wave components are necessarily considered
to increase model accuracy in the multidirectional case as in the
unidirectional case.

Keywords: Ocean waves, Phase-resolved model, Real-time prediction

1. INTRODUCTION

In the study of marine science and ocean engineering, real-
time ocean wave predictions have an impact on operations of
surface vessels [1, 2], ocean wave energy harvesting systems
[3, 4], and designs of offshore structures [5]. In particular, the
control strategies of floating offshore wind turbines (FOWTs)
require an accurate real-time forecast of the incoming surface
wave conditions.

The disadvantage of phase-averaged models (e.g., SWAN,
[6]) in ocean wave prediction is that the lack of phase infor-
mation in the average quantities renders them to have inherent
difficulties in capturing accurate and detailed nonlinear wave dy-
namics. Alternatively, phase-resolved models accounting for a
history of phase correlations provide the instantaneous state of
ocean wave motion and the statistical properties in relation to
surface wave dynamics. In an attempt to forecast the ocean wave
surface from the spatio-temporal, several authors have reported
on the development of phase-resolved reconstruction and predic-
tion algorithms from large spatio-temporal data sets of surface
wave elevations obtained by an optical system [1, 2, 7, 8].

In the present study, we extend the phase-resolved real-time
ocean wave prediction in the case of multidirectional sea states.
Due to the time constraints in the real-time prediction, the wave
model based on the linear wave theory is mainly employed since
the wave prediction model leads to the process of wave model
inversion where we obtain the initial condition as a form of am-
plitude coefficients from the wave elevation measurements in time
and space. However, in real conditions, the inclusion of nonlin-
earity is important to provide a more accurate prediction of wave
surface. The nonlinearity must be included in the simulation of
ocean wave surface, particularly when it comes to the severe sea
states in terms of wave steepness.

In order to achieve a good balance between the computational
efficiency and the model accuracy, we employ wave models on
the basis of Lagrangian analysis, since the Lagrangian model at
the equivalent order, compared to the Eulerian counterparts, ap-
pears to be particularly attractive in the study of waves in steep
conditions [9]. Further, in the case of multidirectional wave
fields, more challenging computational costs are expected due
to the great increases in wave components and wave measure-
ments to address directional spreading in comparison with the
unidirectional case. Hence, we propose a simplified and succinct
assimilation method in the process of reconstruction, which leads
to a significant reduction in the computational effort. The exper-
imental work was carried out at École Centrale de Nantes (ECN) to mimic the measurements of a LIDAR in a real configuration with multidirectional seas (short-crested waves) as well as uni-directional seas (long-crested waves). Through the comparison of the model results by the linear model and the second-order nonlinear model in the unidirectional and multidirectional cases to the experimental data, we investigate the effect of spreading and multidirectional wave effects on real-time wave prediction.

## 2. WAVE MODELS

### 2.1 Linear wave model

A Cartesian coordinate system $(x, y, z) = (\mathbf{r}, z)$ is selected, with $x, y$ located on the mean water surface, and $z$ taken positive vertically upwards from the still water level. For inviscid, incompressible, and irrotational fluid, a linear ocean surface elevation is derived from the linearized water wave boundary problem with respect to wave steepness:

$$ \eta^{\text{LWT}}(x, y, t) = \sum_{n=1}^{N} [a_n \cos \psi_n + b_n \sin \psi_n] \quad (1) $$

where subscript $n$ indicates the $n$th wave component vector, $n = (n_\omega, n_\theta)$ with $n_\omega$ and $n_\theta$ indicating $n_\omega$th components in frequency and direction, respectively. The total number of wave components is $N = N_\omega \times N_\theta$ and a combination of the complex amplitude $A_n$ and the phase $\varphi_n$ yields the wave amplitude parameters $(a_n, b_n) = (A_n \cos \varphi_n, A_n \sin \varphi_n)$. The linear phase function is given:

$$ \psi_n = \mathbf{k}_n \cdot \mathbf{r} - \omega_n t \quad (2) $$

where the propagating direction $\theta_{n\theta}$ with respect to the $+x$-direction determines the unit wave vector wavenumber vector $\hat{\mathbf{k}}_{n\theta} = (\cos \theta_{n\theta}, \sin \theta_{n\theta})$. The wavenumber vector $\mathbf{k}_n$ and the unit wave vector $\hat{\mathbf{k}}_{n\theta}$ are related as $\hat{\mathbf{k}}_{n\theta} = \mathbf{k}_n/k_{n\omega}$. The magnitude of wavenumber $k_{n\omega}$ is related to the wave angular frequency $\omega_{n\omega}$ by the linear dispersion relation in deep water. We note here that the extension to finite water depth is straightforward.

### 2.2 Improved choppy wave model

In order to achieve an efficient simulation, Nouguier et al. [10] derived the choppy wave model (CWM) based on the Lagrangian approach on the nonlinear ocean surface. Since then, to overcome a significant drawback of CWM regarding the lack of nonlinearity in its celerity, Guérin et al. [11] formulated the improved choppy wave model (ICWM) from the second-order Lagrangian expression. To connect the wave information recorded in the form of time series based on the Eulerian coordinate system with the Lagrangian-based models, ICWM was transformed into the approximated model in the Eulerian system [8, 11]. The free surface elevation by ICWM is given as:

$$ \eta^{\text{ICWM}}(x, y, t) = \sum_{n=1}^{N} [a_n \cos \Psi_n + b_n \sin \Psi_n + \frac{1}{2}(a_n^2 + b_n^2)k_n] \quad (3) $$

in which the nonlinear phase function for ICWM $\Psi_n$ retains the nonlinear phase shift and Stokes drift $\mathbf{U}_{s0}$:

$$ \Psi_n = \mathbf{k}_n \cdot [\mathbf{r} - \sum_{i=1}^{N} \mathbf{k}_i (-a_i \sin \tilde{\psi}_i + b_i \cos \tilde{\psi}_i)] - \tilde{\omega}_n t \quad (4) $$

$$ \mathbf{U}_{s0} = \sum_{n=1}^{N} (a_n^2 + b_n^2) \omega_n \mathbf{k}_n \quad (5) $$

where a tilde superscript refers to the modified variables with Stokes drift (i.e., $\tilde{\omega}_n = \omega_n + \frac{1}{2}\mathbf{k}_n \cdot \mathbf{U}_{s0}$ and $\tilde{\psi}_n = \mathbf{k}_n \cdot \mathbf{r} - \tilde{\omega}_n t$). For the zero-mean sea level, the last term in Eq. (3) is introduced.

## 3. DATA ASSIMILATION

Real-time phase-resolved ocean wave predictions are carried out by first data assimilation to reconstruct wave amplitude parameters on the basis of measurements $\tilde{\eta}$, and then wave prediction to simulate wave surfaces over a prediction zone. The wave amplitude parameters $(a_n, b_n)$ are specified by the model inversion from measurements on the assumption of unchanged amplitudes over a certain zone in the space-time domain. By means of a variational approach [12], we determine the optimized parameters via an optimization process minimizing a quadratic cost function representing the difference between predicted and measured ocean surface elevations ($\eta_l$ and $\tilde{\eta}_l$):

$$ F(\mathbf{p}) = \frac{1}{2} \sum_{l=1}^{L} [\eta_l(\mathbf{p}) - \tilde{\eta}_l]^2 = \frac{1}{2} \sum_{j=1}^{J} \sum_{k=1}^{K} [\eta_{(j,k)}(\mathbf{p}) - \tilde{\eta}_{(j,k)}]^2 \quad (6) $$

where the $l$th spatio-temporal measurement is indicated by $l$, a pair of $l = (j, k)$; $j$ and $k$ refer to the measured data at spatial point $r_j = (x_j, y_j)$ and time $t_k$, respectively. We define the total number of spatio-temporal data points $L$ by the size of data in space and time $J$ and $K$, respectively (i.e., $L = J \times K$). Then, the derivatives of the cost function with respect to the parameters $\mathbf{p}$ yield the system of equations with matrix operation $\mathbf{A}\mathbf{p} = \mathbf{B}$:

$$ \begin{aligned} \frac{\partial F}{\partial a_m} &= 0 \Rightarrow \sum_{l=1}^{L} \eta_l \frac{\partial \eta_l}{\partial a_m} = \sum_{l=1}^{L} \tilde{\eta}_l \frac{\partial \eta_l}{\partial a_m} &\Rightarrow A_{(m,n)} p_n &= B_m \\ \frac{\partial F}{\partial b_m} &= 0 \Rightarrow \sum_{l=1}^{L} \eta_l \frac{\partial \eta_l}{\partial b_m} = \sum_{l=1}^{L} \tilde{\eta}_l \frac{\partial \eta_l}{\partial b_m} &\Rightarrow A_{(N+m,n)} p_n &= B_{N+m} \end{aligned} \quad (7) $$

where the model coefficient $1 \times 2N$ vector $\mathbf{p}$ consists of $p_n = a_n$ and $p_{N+n} = b_n$ for $n, m \in \{1, ..., N\}^2$.

### 3.1 Linear assimilation

The assimilation matrices ($\mathbf{A}$ and $\mathbf{B}$) during the wave reconstruction procedure rely on the wave propagating model and are given for LWT:

$$ \begin{aligned} A_{(m,n)}^{\text{LWT}} &= \sum_{l=1}^{L} \cos \psi_{nl} P_{ml}^{\text{LWT}}, & A_{(m,N+n)}^{\text{LWT}} &= \sum_{l=1}^{L} \sin \psi_{nl} P_{ml}^{\text{LWT}} \\ A_{(N+m,n)}^{\text{LWT}} &= \sum_{l=1}^{L} \cos \psi_{nl} Q_{ml}^{\text{LWT}}, & A_{(N+m,N+n)}^{\text{LWT}} &= \sum_{l=1}^{L} \sin \psi_{nl} Q_{ml}^{\text{LWT}} \\ B_m^{\text{LWT}} &= \sum_{l=1}^{L} \tilde{\eta}_l P_{ml}^{\text{LWT}}, & B_{N+m}^{\text{LWT}} &= \sum_{l=1}^{L} \tilde{\eta}_l Q_{ml}^{\text{LWT}} \end{aligned} \quad (8) $$

with

$$ P_{ml}^{\text{LWT}} = \cos \psi_{ml}, \quad Q_{ml}^{\text{LWT}} = \sin \psi_{ml} \quad (9) $$

where $\psi_{ml} = \mathbf{k}_{ml} \cdot \mathbf{r}_l - \omega_m t_l$.
## 3.2 Nonlinear assimilation

Unlike the linear system, since the nonlinear phase function for ICWM $\Psi_n$ containing the nonlinear phase shift and Stokes drift $U_{s0}$ includes the amplitude parameter, we need to proceed iteratively to find the amplitude parameters. The system of nonlinear assimilation was previously derived [8]:

$$ 
\begin{aligned}
A_{(m,n)}^{\text{ICWM}} &= \sum_{l=1}^{L} (\cos \Psi_{nl} + \frac{1}{2} a_n k_n) P_{ml}^{\text{ICWM}} \\
A_{(m,N+n)}^{\text{ICWM}} &= \sum_{l=1}^{L} (\sin \Psi_{nl} + \frac{1}{2} b_n k_n) P_{ml}^{\text{ICWM}} \\
A_{(N+m,n)}^{\text{ICWM}} &= \sum_{l=1}^{L} (\cos \Psi_{nl} + \frac{1}{2} a_n k_n) Q_{ml}^{\text{ICWM}} \\
A_{(N+m,N+n)}^{\text{ICWM}} &= \sum_{l=1}^{L} (\sin \Psi_{nl} + \frac{1}{2} b_n k_n) Q_{ml}^{\text{ICWM}} \\
B_m^{\text{ICWM}} &= \sum_{l=1}^{L} \bar{\eta}_l P_{ml}^{\text{ICWM}}, \quad B_{N+m}^{\text{ICWM}} = \sum_{l=1}^{L} \bar{\eta}_l Q_{ml}^{\text{ICWM}}
\end{aligned}
\quad (10) $$

with

$$ 
\begin{aligned}
P_{ml}^{\text{ICWM}} &= \cos \Psi_{ml} - [k_m(a_m \sin \Psi_{ml} - b_m \cos \Psi_{ml}) \\
&\quad \times \{\sin \tilde{\psi}_{ml} - [k_m(a_m \cos \tilde{\psi}_{ml} + b_m \sin \tilde{\psi}_{ml}) + 1] \\
&\quad \times a_m \omega_m k_m t_l\} + a_m k_m \\
Q_{ml}^{\text{ICWM}} &= \sin \Psi_{ml} - [k_m(a_m \sin \Psi_{ml} - b_m \cos \Psi_{ml}) \\
&\quad \times \{-\cos \tilde{\psi}_{ml} - [k_m(a_m \cos \tilde{\psi}_{ml} + b_m \sin \tilde{\psi}_{ml}) + 1] \\
&\quad \times b_m \omega_m k_m t_l\} + b_m k_m
\end{aligned}
\quad (11) $$

where $\Psi_{ml} = k_{ml} \cdot [\mathbf{r}_l - \sum_{i=1}^{N} (-a_i \sin \tilde{\psi}_{il} + b_i \cos \tilde{\psi}_{il})] - \tilde{\omega}_m t_l$ and $\tilde{\psi}_{il} = k_i \cdot \mathbf{r}_l - \tilde{\omega}_i t_l$.

The solution matrix at the qth iteration is referred to as $\mathbf{p}^q$. The solution at the first iteration is initialized from the linear solution, and the parameter vector $\mathbf{p}^{q+1}$ at the current iteration $q+1$ is obtained based on the computed matrix by using $\mathbf{p}^q$ at the previous iteration $q$. If the relative error between solutions of $\mathbf{p}$ at two successive iterations becomes smaller than a fixed tolerance parameter (set to $10^{-3}$), the assimilation process stops. Also, we consider the maximum iteration number of 100 to account for the severely ill-conditioned or divergent case.

In the matrix system for ICWM shown in Eq. (10), the computational effort is more associated with the second factors (i.e., $P$ and $Q$) rather than the first factors (e.g., $\cos \Psi_{nl} + \frac{1}{2} a_n k_n$). This is due to the components retaining the derivative of the nonlinear phase function with respect to the amplitude coefficients and thus being more variable at every iteration. Accordingly, we can simplify the nonlinear assimilation matrices (A and B) for ICWM by truncating $P$ and $Q$ up to the order of unity. The simplified assimilation nonlinear matrices for ICWM (ICWM-S) consist of

the first parts of ICWM, and $P$ and $Q$ of LWT:

$$ 
\begin{aligned}
A_{(m,n)}^{\text{ICWM-S}} &= \sum_{l=1}^{L} (\cos \Psi_{nl} + \frac{1}{2} a_n k_n) P_{ml}^{\text{LWT}} \\
A_{(m,N+n)}^{\text{ICWM-S}} &= \sum_{l=1}^{L} (\sin \Psi_{nl} + \frac{1}{2} b_n k_n) P_{ml}^{\text{LWT}} \\
A_{(N+m,n)}^{\text{ICWM-S}} &= \sum_{l=1}^{L} (\cos \Psi_{nl} + \frac{1}{2} a_n k_n) Q_{ml}^{\text{LWT}} \\
A_{(N+m,N+n)}^{\text{ICWM-S}} &= \sum_{l=1}^{L} (\sin \Psi_{nl} + \frac{1}{2} b_n k_n) Q_{ml}^{\text{LWT}} \\
B_m^{\text{ICWM-S}} &= \sum_{l=1}^{L} \bar{\eta}_l P_{ml}^{\text{LWT}}, \quad B_{N+m}^{\text{ICWM-S}} = \sum_{l=1}^{L} \bar{\eta}_l Q_{ml}^{\text{LWT}}
\end{aligned}
\quad (12) $$

The simplified method could be justified for use in the procedure of nonlinear assimilation as the derivation of the numerical wave model is based on the assumption of small wave amplitudes. We will verify this simplified approach by comparing the resulting ocean wave surfaces from both assimilation methods (i.e., the previous and simplified ones) with the experimental results later.

## 4. EXPERIMENTAL DATA

We present an experimental setup mimicking measurements coming from remote measurements such as LIDAR (Light detection and ranging) cameras during one of the campaigns of the FLOATECH project. The experiment was conducted at the geometric scale of 1:40 in the hydrodynamic and ocean engineering tank at École Centrale de Nantes (ECN), which is a 30-m wide, 50-m long, and 5-m deep water testing volume. Waves are generated by a wavemaker at one end of the tank, consisting of 48 individual hinged flaps, and absorbed by a beach at the other end. We generated a random wave field using a Pierson-Moskowitz spectrum [13] with a peak period $T_p = 12$ s, a significant wave height $H_s = 7$ m, a peak wavelength $L_p = 225$ m at full scale (corresponding to a moderate nonlinear sea state, $H_s/L_p = 3.1\%$ and deep water, $k_p d \approx 5.6$), and its main propagating direction along the x-axis, $0^\circ$. The wave directionality is considered by the angular spreading function proposed by Mitsuyasu et al. [14]:

$$ G(\theta) = \frac{2^{2s-1} (s!)^2}{180 (2s)!} \cos^{2s} \left(\frac{\theta}{2}\right) \quad \text{for } \theta \in [-180^\circ, 180^\circ] \quad (13) $$

Figure 1 details an experimental setup for a multidirectional wave field obtained from the optical remote sensing system where the sensor is located at a location of $(x_c, y_c, z_c) = (19.86$ m, 15 m, 30 m) and aimed at some distance ($\alpha = 76^\circ$) with vertical and horizontal aperture angles ($\alpha_a = 20^\circ$ and $\beta_a = 64^\circ$). Through $J_r \times J_\theta$ rays with $J_r = 20$ and $J_\theta = 9$, the measurement zone forms a shape of a radial segment (see Figure 2). The wave measurements at all the wave gauges in the experiment satisfied repeatability. We reconstructed a network of wave gauges thanks to the repeatability as well as the rotating straight structure. We made this choice of probe locations in this experiment without considering the actual locations that would result from the intersection between optical rays and an irregular ocean surface.
<figure><img src="image_3.png" /><figcaption>FIGURE 1: EXPERIMENTAL SETUP BY AN OPTICAL SYSTEM (●: WAVE OBSERVATIONS; ▲: TARGET LOCATION, WG21; ■: THREE ADDITIONAL DOWNSTREAM WG22, WG23, AND WG24).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2: LOCATION OF WAVE GAUGES (●: WAVE OBSERVATIONS; ▲: TARGET LOCATION, WG21; ■: THREE ADDITIONAL DOWNSTREAM WG22, WG23, AND WG24).</figcaption></figure>

Of several wave conditions in the campaign of the FLOAT-ECH project, we address case A with moderate directional spreading and wave steepness ($s = 25$ and $H_s/L_p = 3.1\%$). Figure 3 shows the normalized spectral energy density of the free surface elevation $S_\eta^* = S_\eta f_p / (H_s^2/16)$. The straight ladder on which 20 wave gauges are located, with the gaps between gauges varying with the distance from the center of rotation, was allowed to rotate with the evenly-spaced nine angles ($J_\theta = 9$) on the xy-plane ranging from -32° to 32°. Also, as a reference case, the second dataset used to analyze the performance of wave models is a unidirectional case B with the same values of $H_s$ and $T_p$ where a structure parallel to the x-axis was used only and $s = \infty$ denotes the long-crested waves in Table 1. Four additional downstream wave gauges were installed along $y = y_c$, including the target location (referred to as WG21) and farthest gauge WG24 for the

gauge 0.12$L_p$ far away from WG21.

<figure><img src="image_5.png" /><figcaption>FIGURE 3: OCEAN WAVE DIRECTIONAL SPECTRUM OF CASE A.</figcaption></figure>

<table><caption>TABLE 1: EXPERIMENTAL SETUP AND WAVE PARAMETERS.</caption><thead><tr><th>Case</th><th>s</th><th>J<sub>θ</sub></th><th>H<sub>s</sub> (m)</th><th>T<sub>p</sub> (s)</th></tr></thead><tbody><tr><td>A</td><td>25</td><td>9</td><td>7</td><td>12</td></tr><tr><td>B</td><td>∞</td><td>1</td><td>7</td><td>12</td></tr></tbody></table>

# 5. PREDICTION ZONE

## 5.1 Cutoff frequencies and directions

To use the dataset obtained from the experimental test campaign, a judicious choice of cutoff frequencies and directions should be made in advance to provide a suitable phase-resolved wave prediction. As in the previous study [8], we determine the cutoff frequencies based on the spectral density at peak angular frequency $\omega_p$ with a small parameter of $\mu = 0.05$:

$$S_{\eta}(\omega_{\min}) = S_{\eta}(\omega_{\max}) = \mu S_{\eta}(\omega_{p}) \quad (14)$$

where $S_\eta(\omega)$ refers to the wave spectrum at the angular frequency $\omega$, and $\omega_{\min}$ and $\omega_{\max}$ yields the group velocities $c_{gmax}$ and $c_{gmin}$ by the deep water approximation of linear dispersion relation, respectively. Next, in order to retain most of the spectral energy, we set the directional cutoff limits $\theta_{\min} = -45^\circ$ and $\theta_{\max} = 45^\circ$; this directional range contains 99.5% of the spectral energy in the directional case A (see Figure 3).

## 5.2 Temporal evolution of prediction zone

To address the prediction zone over which the wave prediction is available, a similar approach based on the group velocity within the prespecified frequency and direction bandwidths is
followed as in the previous study [8, 15, 16]. We here suggest the prediction zone when measurements are provided using an optical system mounted on the offshore structure over an assimilation time $T_a$. The multidirectional predictable spatio-temporal domain is determined by the intersection of a unidirectional prediction zone within the prescribed wave propagating direction bandwidth $\theta \in [-45^\circ, 45^\circ]$.

Figure 4 illustrates the observation zone and the prediction zone at the latest time of wave reconstruction $t = t_r$. For the sake of simplicity, we consider the simplified boundaries of observation by using the left and right boundaries with the selected measurement points.

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 4: OBSERVATION AND PREDICTION ZONES: (A) OBSERVATION ZONE; (B) PREDICTION ZONE (DASHED: LEFT BOUNDARY; SOLID: RIGHT BOUNDARY; RED SHADING: NOWCAST; GRAY SHADING: FORECAST).</figcaption></figure>

The segment of each boundary propagates in +x-direction differently depending on its propagating direction $\theta$, therefore, the greater the propagating direction $\theta$, the faster the segment propagates in +x-direction. As a result, the temporal prediction zone $t' = t - t_r \in [t'_{min}, t'_{max}]$ at any location $x$ along $y = y_c$ is

given as:

$$ 
\begin{aligned}
t'_{\min} &= \max[-T_a, -T_a + \frac{x - x_H}{c_g \min}] \\
t'_{\max} &= \min[\frac{x - x_D}{c_g \max} \cos(\theta_B/2), \frac{x - x_B + y_B}{c_g \max} \cos 45^\circ]
\end{aligned} 
\quad (15) 
$$

where the half of horizontal aperture angle is $\theta_B = \beta_a/2$. Further, by substituting $\theta_B = 0$, the temporal prediction zone for the multidirectional case is reduced to that of the unidirectional case:

$$ \max[-T_a, -T_a + \frac{x - x_H}{c_g \min}] \le t' \le \frac{x - x_D}{c_g \max}. \quad (16) $$

# 6. NUMERICAL SIMULATIONS

## 6.1 Normalized misfit error

In order to provide an accurate evaluation of model performance, we employ an ensemble average of normalized misfit error by using partly overlapping surface samples, shifted in time by $\Delta t$:

$$ \varepsilon(x, y, t) = \frac{1}{N_s} \sum_{i=1}^{N_s} \frac{|\eta_{pred,i}(x, y, t) - \eta_{ref,i}(x, y, t)|}{H_s} \quad (17) $$

where $\eta_{pred,i}$ and $\eta_{ref,i}$ denote the prediction and reference values of surface elevation from the ith surface sample, respectively.

Here, the time shift between two consecutive surface samples $\Delta t$ plays a significant role in the total number of surface samples $N_s$ yielding the convergence. We use the number of surface samples $N_s = 400$ (or $T_c/T_p \approx 25$) which is required to reach the convergence in prediction error for both cases. The smaller $\Delta t$, the faster convergence for the prediction error, resulting in the less total time duration during the assimilation procedure over $N_s$ samples for the convergence (i.e., $T_c = T_a + (N_s - 1)\Delta t$) despite the convergence error which is independent of $\Delta t$ [8]. Accordingly, we use $\Delta t/T_p = 0.053$ which is small enough to ensure converged results. The normalized misfit error is further averaged over the time prediction zone $[t_{min}, t_{max}]$:

$$ \varepsilon^p(x, y) = \frac{1}{t_{\max} - t_{\min}} \int_{t_{\min}}^{t_{\max}} \varepsilon(x, y, t) dt \quad (18) $$

## 6.2 Wave parameters

In the following, we employ the spatio-temporal ocean surface (where the number of data is $J \times K$) acquired by an optical sensor at grazing incidence in the process of wave assimilation/reconstruction, and then evolve in time and space to offer the ocean surface forecast. We used all the measured data for each experimental case, therefore, the numbers of measurement points in space $J = J_r \times J_\theta$ for Cases A and B are 20×9 and 20, respectively. The assimilation time $T_a$, over which the wave field is reconstructed from the wave measurement, directly yields the number of measurement data in time $K$. In order to gain an efficient numerical effort and model accuracy simultaneously, we conducted an optimization process with respect to the assimilation time $T_a$ as shown in Figure 5. For all the cases of this study, the normalized misfit errors are converged for $T_a/T_p \approx 5.2$, corresponding to $K = 100$.
<figure><img src="image_8.png" /><figcaption>FIGURE 5: PREDICTION ERROR OF ICWM AT WG21 AGAINST $T_a/T_p$ (SOLID: CASE A; DASHED: CASE B).</figcaption></figure>

The bandwidths of the reconstructed wave field are required to be defined from the numerical point of view, where the wavenumber bandwidths are determined by the advection of wave information during $T_a$ as well as the distance between wave gauges (or observation grid). Accordingly, the wavenumber of the first frequency mode is obtained by the sum of the largest horizontal distance with respect to the x-axis between wave gauges, namely, the beginning and end of gauges ($x_b$ and $x_e$) and the propagating distance of the right boundary over $T_a$ into +x-direction (see Figure 4):

$$k_1 = \frac{2\pi}{x_e - x_b + c_{g,N_\omega} T_a} \quad (19)$$

Here, it confirms that the last wave component in frequency ($n_\omega = N_\omega$) is necessarily prespecified in defining the minimum cutoff limit (or $k_1$). Following Desmars et al. [8], we use a linear proportional relation between $k_{N_\omega}$ and $N_\omega$ with $k_{N_\omega} = 20k_p$ and $N_\omega = 50$ in order to fix the identical spacing in between the cutoff limits. The prediction error decreases with increasing $N_\omega$, and converges to a minimum at $N_\omega = 30$, corresponding to $k_{N_\omega} = 12k_p$, which is consistent with the finding of Desmars et al. [7] that $k_{N_\omega} = 12k_p$ may be sufficiently large for addressing an appropriate inversion of the wave reconstruction. For the directional bandwidth, we use the same limits as in the prediction zone [$\theta_{min} = -45^\circ, \theta_{max} = 45^\circ$] in the numerical simulation. As for the number of wave components in frequency $N_\omega$, the sensitivity analysis with respect to the number of directional wave components $N_\theta$ is also conducted.

Figure 6 shows that the optimal number of wave components for multidirectional case A is determined as $N = N_\omega \times N_\theta = 30 \times 7$ from the sensitivity analysis. The dash-dot line denotes the results for the unidirectional case B by using the identical number of wave components in frequency, that is, $N = N_\omega \times N_\theta = 30$, which can provide a precise comparison with a consistent approach in the unidirectional case. While the change in the frequency bandwidth with respect to $N_\omega$ makes little difference in results, the increase of $N_\theta$ leads to a sharp reduction in the prediction error in comparison to the unidirectional approach by $N_\theta = 1$. It is confirmed from the comparison of results between both cases that the multidirectional approach including the directional aspect reaches the model accuracy which is equivalent to that in the unidirectional case.

<figure><img src="image_9.png" /><figcaption>FIGURE 6: PREDICTION ERROR OF ICWM AT WG21: (A) AGAINST $N/\theta$; (B) AGAINST $N/\omega$ (SOLID: CASE A; DASHED: CASE B).</figcaption></figure>

## 6.3 Results

First, the ability of the simplified nonlinear assimilation is verified by comparing the prediction errors of the previous and the newly-developed methods. The time length of the practical prediction zone from the latest time of nowcast $t_r$ to when the reconstructed information becomes obsolete $t_{max}$ at WG21, namely, $t'_{max} = t_{max} - t_r$ is about $1.8T_p = 21.6$ s at full scale. Therefore, the constraint for a real-time system is that the total numerical time $T_{total}$, summation of times for nowcast $T_{nowcast}$ and forecast $T_{forecast}$, has to be smaller than $1.8T_p = 21.6$ s. Table 2 compares the prediction error ($\epsilon^p$) as well as numerical cost and efficiency. When it comes to the time duration, the nowcast reconstructing the wave amplitudes far outweighs the forecast simulating wave surfaces with the reconstructed wave amplitudes (i.e., $T_{nowcast} \gg T_{forecast}$). The simplified method appears to offer a significant improvement with respect to computational efficiency and stability over the previous one, with maintaining the model performance. In particular, the simplified method allows real-time phase-resolved wave forecasting ($T_{total} < t'_{max} = t_{max} - t_r$). Here, it would be opportune to discuss the benefit of reducing the computational cost. The final objective of the wave prediction is to be able to control the wind turbine to minimize the motions of the platform and/or to maximize energy production. The wind turbine controller needs to know the incoming waves or the resulting excitation force in advance in order to perform the control. With a lead time of 16.2 s, this should allow the set-up of such a control strategy on a FOWT. We notice that this significant improvement is consistent with the increase in the spatial extent of the prediction zone that can have some practical applications.

Figure 7 compares the descriptions of wave surface by ICWM and LWT with laboratory observations. Model predictions are
<table><caption>TABLE 2: COMPARISON BETWEEN ASSIMILATION METHODS FOR CASE A AT WG21 WHERE PRACTICAL PREDICTION ZONE IN TIME $t_{max} - t_{r} \approx 21.6$ S AT FULL SCALE.</caption><thead><tr><th>Method</th><th>Previous, Eq.(10)</th><th>Simplified, Eq.(12)</th></tr></thead><tbody><tr><td>$T_{nowcast}$</td><td>59.04</td><td>5.40</td></tr><tr><td>$T_{forecast}$</td><td>&lt; 0.01</td><td>&lt; 0.01</td></tr><tr><td>$T_{total}$</td><td>59.04</td><td>5.40</td></tr><tr><td>Iteration number (max = 100)</td><td>23.3</td><td>6.2</td></tr><tr><td>Convergence probability ($N_s = 400$)</td><td>88%<br>(346)</td><td>97%<br>(388)</td></tr><tr><td>$\epsilon^p$</td><td>0.0429</td><td>0.0431</td></tr></tbody></table>

in overall good agreement with the observations, both for multidirectional and unidirectional cases. In order to quantitatively assess the model performances in the multidirectional case A, both models are further checked with the temporal evolution of the ensemble average of misfit error $\epsilon$ over the surface samples within the temporal prediction zone [$t_{min}, t_{max}$] (see Figure 8). As shown in Figure 8(a), the theoretical prediction zone overlaps with the time horizon where the prediction errors remain significantly small, where the value of minimum error is 3.99% by LWT, which is decreased by about 9% when using ICWM. Moreover, ICWM is somewhat more advantageous at $t = t_{max}$ with the maximum misfit error of 5.93%, while it is 6.74% by LWT, which is about 14% larger compared to ICWM. Further, Figure 8(b) shows that the theoretical prediction zone in time and space agrees successfully with the resulting surface error at all the locations along the centerline.

<figure><img src="image_10.png" /><figcaption>FIGURE 7: TIME SERIES OF SURFACE ELEVATION AT WG21: (A) CASE A; (B) CASE B (SOLID: DATA; DASHED: ICWM; DOTTED: LWT).</figcaption></figure>

Lastly, we examine the effect of the directional spreading and the propagating distance from the wave measurement zone on the phase-resolved wave forecast (see Figure 9). Above all, the model performances by both linear and nonlinear models become relatively more deleterious with the increasing distance from the observation zone. In other words, the most significant prediction errors in both cases become the greatest for all the models at WG24, to which the propagating distances from the end of $x_e$ are about $0.28L_p$. It is not a surprising result since the wave

<figure><img src="image_11.png" /><figcaption>FIGURE 8: EVOLUTION OF PREDICTION ERROR FOR CASE A: (A) USING ICWM AND LWT AT WG21; (B) USING ICWM AT ALL WGS ALONG $y = y_c$ (SOLID: ICWM; DASHED: LWT; DASH-DOT: BOUNDARIES OF PREDICTION ZONE; RED SHADING: NOWCAST; GRAY SHADING: FORECAST; VERTICAL WHITE LINES: x-LOCATION OF WGS ALONG $y = y_c$; BLACK RECTANGLE: ASSIMILATED DATA SET).</figcaption></figure>

information obtained from the system of equations to assimilate observation data is increasingly in error as waves propagate in space and time.

Regardless of whether the sea state is unidirectional or directional, the nonlinear models demonstrate nearly the same results at every downstream location in terms of the normalized misfit error. In contrast, the ratio of prediction errors by ICWM to that by LWT is quite different between the multidirectional and unidirectional cases; the performance improvement by the nonlinear model is more pronounced for the unidirectional cases. In comparison to the model based on linear wave theory, the nonlinear model on the basis of Lagrangian analysis retains the summation in Eq. (4) (or geometric nonlinearity) and Stokes drift $U_{s0}$ in its equation. The angular directionality leads the effect of Stokes drift to be negligible in comparison to the unidirectional case. Stokes drift in the nonlinear phase function is represented as a
second-order term in wave steepness, therefore, the summation of the directional wave components becomes less dominant than that in the unidirectional condition. This is mainly due to the fact that the summation of the more decomposed components in the direction for one frequency mode is smaller.

<figure><img src="image_12.png" /><figcaption>FIGURE 9: SPATIAL EVOLUTION OF PREDICTION ERROR: (A) ERROR BY ICWM; (B) COMPARISON OF ERROR BETWEEN ICWM AND LWT (SOLID: CASE A; DASHED: CASE B).</figcaption></figure>

## 7. CONCLUSION

The wave prediction algorithm with the simplified assimilation method is validated against the dedicated experimental data in the directional case with maintaining the model accuracy as in the unidirectional case. The major findings are given below: (1) the advantage of using the simplified assimilation approach is the significant reduction in numerical cost; (2) the overall results demonstrate that the model predictions by LWT and ICWM are in overall good agreement with the observations; (3) the performance improvement by the nonlinear model is less obvious for the directional case due to the decomposition into the directional wave components.

## ACKNOWLEDGMENTS

This work was done within the framework of the European H2020 FLOATECH project, Grant agreement ID: 101007142. The authors acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-20-CE05-0039 (project CREATIF).

## REFERENCES

[1] Grilli, Stéphan T, Guérin, Charles-Antoine and Goldstein, Bart. "Oceanwave reconstruction algorithms based on spatio-temporal data acquired by a flash LiDAR camera."

The Twenty-first International Offshore and Polar Engineering Conference: pp. 275-282. 2011. OnePetro.

[2] Nouguier, Frédéric, Grilli, Stéphan T and Guérin, Charles-Antoine. "Nonlinear ocean wave reconstruction algorithms based on simulated spatiotemporal data acquired by a flash LiDAR camera." *IEEE Transactions on Geoscience and Remote Sensing* Vol. 52 No. 3 (2013): pp. 1761-1771.

[3] Li, Guang, Weiss, George, Mueller, Markus, Townley, Stuart and Belmont, Mike R. "Wave energy converter control by wave prediction and dynamic programming." *Renewable Energy* Vol. 48 (2012): pp. 392-403.

[4] Previsic, Mirko, Karthikeyan, Anantha and Lyzenga, David. "In-ocean validation of a deterministic sea wave prediction (DSWP) system leveraging X-band radar to enable optimal control in wave energy conversion systems." *Applied Ocean Research* Vol. 114 (2021): p. 102784.

[5] Ma, Yu, Sclavounos, Paul D, Cross-Whiter, John and Arora, Dhiraj. "Wave forecast and its application to the optimal control of offshore floating wind turbine for load mitigation." *Renewable Energy* Vol. 128 (2018): pp. 163-176.

[6] Booij, NRRC, Ris, Roeland C and Holthuijsen, Leo H. "A third-generation wave model for coastal regions: 1. Model description and validation." *Journal of geophysical research: Oceans* Vol. 104 No. C4 (1999): pp. 7649-7666.

[7] Desmars, Nicolas, Pérignon, Yves, Ducrozet, Guillaume, Guérin, Charles-Antoine, Grilli, Stephan T and Ferrant, Pierre. "Phase-resolved reconstruction algorithm and deterministic prediction of nonlinear ocean waves from spatiotemporal optical measurements." *International Conference on Offshore Mechanics and Arctic Engineering*, Vol. 51272: p. V07BT06A054. 2018. American Society of Mechanical Engineers.

[8] Desmars, Nicolas, Bonnefoy, Félicien, Grilli, ST, Ducrozet, Guillaume, Perignon, Yves, Guérin, C-A and Ferrant, Pierre. "Experimental and numerical assessment of deterministic nonlinear ocean waves prediction algorithms using non-uniformly sampled wave gauges." *Ocean Engineering* Vol. 212 (2020): p. 107659.

[9] Pierson Jr, Willard J. "Models of random seas based on the Lagrangian equations of motion." Technical report no. New York Univ Bronx School of Engineering and Science. 1961.

[10] Nouguier, Frédéric, Guérin, Charles-Antoine and Chapron, Bertrand. ""Choppy wave" model for nonlinear gravity waves." *Journal of geophysical research: oceans* Vol. 114 No. C9 (2009): pp. 1-16.

[11] Guérin, Charles-Antoine, Desmars, Nicolas, Grilli, Stéphan T, Ducrozet, Guillaume, Perignon, Yves and Ferrant, Pierre. "An improved Lagrangian model for the time evolution of nonlinear surface waves." *Journal of Fluid Mechanics* Vol. 876 (2019): pp. 527-552.

[12] Blondel, E, Bonnefoy, F and Ferrant, P. "Deterministic non-linear wave prediction using probe data." *Ocean Engineering* Vol. 37 No. 10 (2010): pp. 913-926.

[13] Pierson Jr, Willard J and Moskowitz, Lionel. "A proposed spectral form for fully developed wind seas based on the
similarity theory of SA Kitaigorodskii.” *Journal of geo-physical research* Vol. 69 No. 24 (1964): pp. 5181–5190.

[14] Mitsuyasu, Hisashi, Tasai, Fukuzo, Suhara, Toshiko, Mizuno, Shinjiro, Ohkusu, Makoto, Honda, Tadao and Riki-ishi, Kunio. “Observations of the directional spectrum of ocean WavesUsing a cloverleaf buoy.” *Journal of Physical Oceanography* Vol. 5 No. 4 (1975): pp. 750–760.

[15] Wu, Guangyu. "Direct simulation and deterministic prediction of large-scale nonlinear ocean wave-field." Ph.D. Thesis, Massachusetts Institute of Technology. 2004.

[16] Qi, Yusheng, Wu, Guangyu, Liu, Yuming and Yue, Dick KP. "Predictable zone for phase-resolved reconstruction and forecast of irregular waves." *Wave Motion* Vol. 77 (2018): pp. 195–213.