# Nonlinear deterministic reconstruction and prediction of remotely measured ocean surface waves

Nicolas Desmars, Moritz Hartmann, Jasper Behrendt, Norbert Hoffmann, Marco Klein

## ► To cite this version:

Nicolas Desmars, Moritz Hartmann, Jasper Behrendt, Norbert Hoffmann, Marco Klein. Nonlinear deterministic reconstruction and prediction of remotely measured ocean surface waves. Journal of Fluid Mechanics, 2023, 975, pp.A8. 10.1017/jfm.2023.841. hal-04279597

HAL Id: hal-04279597

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

Submitted on 10 Nov 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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>

# Nonlinear deterministic reconstruction and prediction of remotely measured ocean surface waves

Nicolas Desmars<fnref n="1" />,<fnref n="2" />,†, Moritz Hartmann<fnref n="2" />,<fnref n="3" />, Jasper Behrendt<fnref n="4" />, Norbert Hoffmann<fnref n="1" />,<fnref n="5" /> and Marco Klein<fnref n="1" />,<fnref n="2" />

(Received 6 February 2023; revised 27 September 2023; accepted 30 September 2023)

Algorithms for reconstructing and predicting nonlinear ocean wave fields from remote measurements are presented. Three types of synthetic observations are used to quantify the influence of remote measurement modulation mechanisms on the algorithms' performance. First, the observations correspond to randomly distributed surface elevations. Then, they are related to a marine radar model – the second type takes the wave shadowing modulation into account whereas the third one also includes the tilt modulation. The observations are numerically generated based on unidirectional waves of various steepness values. Linear and weakly nonlinear prediction algorithms based on analytical models are considered, as well as a highly nonlinear algorithm relying on the high-order spectral (HOS) method. Reconstructing surfaces from shadowed observations is found to have an impact limited to the non-visible regions, while tilt modulation affects the reconstruction more generally due to the indirect, more complex extraction of wave information. It is shown that the accuracy of the surface reconstruction mainly depends on the correct modelling of the wave shape nonlinearities. Modelling the nonlinear correction of the dispersion relation, in particular the frequency-dependent wave phase effects in the case of irregular waves, substantially improves the prediction. The suitability of the algorithms for severe wave conditions in finite depth and using non-perfect observations is assessed through wave tank experiments. It shows that only the third-order HOS solution predicts the right amplitude and phase of an emerging extreme wave, emphasizing the relevance of the corresponding physical modelling.

† Email address for correspondence: nicolas.desmars@protonmail.com

<footnote n="1">¹Dynamics Group, Hamburg University of Technology, 21073 Hamburg, Germany</footnote>

<footnote n="2">²Ship Performance Department, Institute of Maritime Energy Systems, German Aerospace Center (DLR), 21502 Geesthacht, Germany</footnote>

<footnote n="3">³Institute for Ship Structural Design and Analysis, Hamburg University of Technology, 21073 Hamburg, Germany</footnote>

<footnote n="4">⁴BHS Hamburg, 21073 Hamburg, Germany</footnote>

<footnote n="5">⁵Mechanical Engineering, Imperial College London, SW7 2AZ London, UK</footnote>
# 1. Introduction

The deterministic (i.e. phase-resolved) prediction of the incoming wave loads on marine structures represents one of the main challenges towards more efficient and safer offshore operations. The cost reduction of offshore wind through the extension of operational time windows for the installation/maintenance of offshore wind turbines, the gain in efficiency of marine energy devices through active control or the safety improvement of aircraft takeoffs/landings are only a few examples of domains that could benefit from the advance knowledge of wave-induced structure motions. Such predictions are performed by means of the calculation of the ocean surface evolution from the analysis of surrounding wave measurements. More specifically, they consist of three steps: (i) measuring some wave-related quantities, (ii) reconstructing the underlying wave field (i.e. extracting wave information from the measurements) and (iii) propagating the reconstructed wave field to the area of interest.

Regardless of their well-characterized properties and proven suitability for deterministic predictions, in *situ* wave measurements, such as made from wave buoys (e.g. Fisher, Thomson & Schwendeman 2021) or acoustic Doppler current profilers (Huchet et al. 2021), are not adaptable to predictions over a region in motion, such as around a ship with forward speed, due to the need to update the measured region according to the trajectory of the area of interest. In that case, measurements are typically performed from remote sensors mounted on a moving structure. Lidar cameras, for example, have recently been developed for the measurement of ocean surface displacements at a distance (Belmont et al. 2007; Kabel, Georgakis & Rod Zeeberg 2019), providing suitable data for wave prediction (Grilli, Guérin & Goldstein 2011; Nouguier, Grilli & Guérin 2014). Stereovideo systems are also able, through the mapping of the surface elevation, to give access to sufficient information for short-term forecasts (Mérigaud & Tona 2022). However, the real-time accessibility of wave data from stereo imaging is still constrained by the ‘significant computational time required to extract the three-dimensional elevation maps from a pair of images’ (Guimarães et al. 2020), and subject to strong visibility limitations. Moreover, both lidar cameras and stereovideo systems are limited to measurement distances of ~100 m, restraining the prediction horizon to a few characteristic wave periods. Even if such prediction horizons are already valuable for applications such as optimal control of marine energy harvesting systems (e.g. Li et al. 2012; Ma et al. 2018), using sensors of higher measurement range could benefit a much larger variety of offshore operations.

The most commonly used remote sensors for the analysis of the ocean surface surrounding an offshore structure are non-coherent X-band marine radars, which are capable of measuring ocean surfaces over a much larger spatial domain, that is within a radius of ~3 km around the antenna. Surface wave identification from non-coherent radar data relies on the interpretation of the backscattered intensity which depends, through Bragg resonance, on the slope modulations of short ripples (gravity-capillary waves) by longer gravity waves (e.g. Alpers & Hasselmann 1978; Alpers, Ross & Rufenach 1981). In addition, both the radars and the aforementioned remote sensing instruments located onboard a structure acquire data at grazing incidence. They are thus subject to wave shadowing that creates spatial gaps in the measured data sets. Two main approaches are used to retrieve the wave field from radar measurements (see the extensive review of
Nonlinear prediction of remotely measured ocean waves

Huang, Liu & Gill 2017). The traditional method employs a semi-empirical modulation transfer function that relates the modal amplitudes of radar intensity to those of surface elevation in the frequency–wavenumber space (Young, Rosenthal & Ziemer 1985; Nieto Borge et al. 2004). Because this method identifies gravity wave components by choosing those which satisfy the linear dispersion relation, nonlinearities are disregarded. The alternative approach relies on the direct resolution of the modulation mechanisms of the radar backscattered intensity by assuming its proportionality to the surface slope. In that case, the linear wave theory is also assumed to facilitate data processing (Dankert & Rosenthal 2004; Naaijen et al. 2018; Simpson et al. 2020).

Nevertheless, the linear wave assumptions have been shown to limit the prediction accuracy as the wave steepness and prediction horizon increase. To overcome this issue, second-order wave models that better represent the steep waves' sharpness have been investigated (e.g. Hlophe et al. 2021), albeit with a limited impact on the prediction accuracy because they miss nonlinear phase effects. Including corrections of the dispersion relation, appearing at third order, was shown to significantly help improve predictions (e.g. Blondel-Couprie, Bonnefoy & Ferrant 2013; Desmars et al. 2020; Meisner et al. 2023), as well as using higher-order models based on the modified nonlinear Schrödinger equation (e.g. Trulsen 2007; Simanesew et al. 2017). For large wave steepness, the high-order spectral (HOS) method theoretically outperforms the aforementioned wave models (Klein et al. 2020) by considering the evolution of the surface dynamics up to any arbitrary order of nonlinearity, with a relatively high computational efficiency. However, despite the rather large number of studies on HOS predictions, either based on sequential (Yoon, Kim & Choi 2016; Wang & Pan 2021; Wang et al. 2022) or variational (Aragh & Nwogu 2008; Blondel, Bonnefoy & Ferrant 2010; Qi et al. 2018a; Fujimoto & Waseda 2020; Wu, Hao & Shen 2022) assimilation strategies, the compatibility of the developed algorithms with remote measurement techniques must be further evaluated. Indeed, previous works suppose that the measurements (i) provide direct access to the surface elevation and (ii) are uniformly sampled in space and time. These two assumptions are not verified by non-coherent X-band radars in realistic conditions, for which the measurements take the form of return intensities subject to (a priori unknown) wave shadowing. Moreover, rather than a generally assumed sequence of spatial snapshots, radar data describe a space/time helix (Al-Ani et al. 2019). The available data are thus indirectly related to surface elevation and distributed in space and time in a non-trivial manner. It is expected that both of these features have an influence on the capacity of the methods to access information about the surface dynamics and predict the wave field.

It is thus proposed in the present paper to assess the impact of the modulation mechanisms that non-coherent radar measurements are subject to on the performance of linear, weakly nonlinear and highly nonlinear prediction algorithms. With this aim, three types of measurements (also called observations) with different levels of modulation are considered as input for the prediction of unidirectional wave fields with various wave steepness values. First, the observations correspond to surface elevations randomly distributed in space and time. Then, they are related to a radar model – the second type takes the shadowing modulation into account whereas the third one also includes the tilt modulation that makes the observed quantity mainly dependent on the surface slope. Regarding the prediction algorithms, two approaches are employed. The first one, based on the inversion of analytical wave models, extends the linear and weakly nonlinear approaches proposed by Desmars et al. (2020) to radar observations and to the assessment of the reconstructed/predicted surface potential. The second one relies on the HOS method and extends the numerical basis of Desmars et al. (2022)
## N. Desmars and others

to irregular waves and to radar observations. It follows a rather classical iterative optimization procedure that seeks the optimal variables (here, the HOS model parameters) by fitting the observations during a predefined time interval, such as done by variational assimilation, but presents some interesting numerical properties that are detailed later in the paper. By inter-comparison of the performance pertaining to each algorithm, the important hydrodynamic properties are highlighted for either the wave field reconstruction or prediction, which are independently evaluated. The case of regular waves is first considered to clearly identify the characteristics of each modelling method, then irregular waves are investigated. High-fidelity synthetic wave data are first used, before applying the methods to experimental wave tank data.

The paper is structured as follows. The theoretical aspects of non-coherent radar measurements, followed by those of the algorithms used for the wave field reconstructions, are presented in § 2. Section 3 details the numerical set-up that is used to evaluate the performance of the surface reconstructions and discusses the results. Likewise, § 4 describes the chosen configuration to study the methods’ prediction performance, before analysing the results. In § 5, the reconstruction/prediction methods are evaluated against wave tank experiments. Finally, § 6 presents the overall conclusions.

## 2. Theoretical background

Predicting surface waves deterministically implies going through three major steps. The first one is the measurement of the quantities that contain relevant information about the wave field for its prediction over a specific area of interest. The second step, referred to as surface reconstruction, is the extraction of wave information from the measurements to initialize the physical model that is used for the propagation of the wave field. The last step is the propagation of the reconstructed surface, giving access to the prediction. Considering that the ocean surface has to be predicted around a potentially mobile structure (e.g. a ship with forward speed), measurements are assumed to be done at a distance by a sensor mounted onboard the structure. The particularities of such remote measurements are detailed in this section, emphasizing their spatial non-uniformity and indirect relationship with surface fields. Then the investigated linear, weakly nonlinear and highly nonlinear methods for surface reconstruction are presented. For the sake of completeness, the theoretical aspects of the methods detailed in this paper describe the general case of directional wave fields, although the applications presented later are restricted to unidirectional waves.

### 2.1. Remote measurements

The prediction of the wave field around a mobile structure has to rely on measurements whose locations adapt to the structure’s trajectory. X-band marine radars are considered in the subsequent developments. In contrast to coherent radars that can provide an estimation of the waves’ orbital velocities through the measurement of the Doppler frequency of the backscattered signal (e.g. Lyzenga et al. 2015; Støle-Hentschel et al. 2018), the considered sensor in the following applications is the much more common non-coherent radar, i.e. it only gives access to the backscattered signal intensity.

#### 2.1.1. Marine radar model

Raw marine radar data are given in a polar coordinate system (r, β, z). As depicted in figure 1, r denotes the range distance, i.e. the distance between the antenna and the
<figure><img src="image_3.png" /><figcaption>Figure 1. Geometrical parameters used for the representation of the shadowing and tilt modulations of marine radar beams.</figcaption></figure>

horizontal location of the ocean surface whose elevation is referred to as $\eta$, and $\beta$ is the azimuth angle of the radar look direction.

The r-axis is located at the mean surface level and points towards the radar look direction, while the z-axis is vertical and positive upward. The range resolution is constant and decided by the radar beam's pulse period, and the azimuth resolution is determined by the antenna beamwidth. For convenience, the ocean surface dynamics is described in a Cartesian system $(x, y, z) = (x_a + r \cos \beta, y_a + r \sin \beta, z)$, with $(x_a, y_a)$ the horizontal coordinates of the radar antenna. Even though many studies consider radar data to be a sequence of instantaneous two-dimensional spatial snapshots (usually called 'radar images'), the spatio-temporal distribution of the radar data actually describes a helix (Al-Ani et al. 2019). For unidirectional waves, however, the helix structure of the data does not apply, and they are modelled as a sequence of one-dimensional snapshots separated by a time step that corresponds to the antenna rotation period.

Among the different modulation mechanisms of the backscattered signal by the ocean surface, the shadowing and tilt modulations are considered prominent (e.g. Naaijen & Wijaya 2014; Salcedo-Sanz et al. 2015) and their respective impact on wave prediction is studied in the present paper. The shadowing modulation describes the impact of obstructing waves in the region illuminated by the radar. Although the actual wave shadowing follows complex diffraction processes (Plant & Farquharson 2012), it can be conveniently approximated by a geometrical approach (see figure 1): a point on the ocean surface is shadowed if a continuous straight line between that point and the radar antenna cannot be drawn without crossing the ocean surface. Practically, the shadowing modulation produces gaps in the spatial distribution of the radar measurements, as shown in figure 2(a). Because shadowing is stronger for larger angles of incidence $\Theta$ (defined in figure 1), the number and size of the gaps increase with the measurement distance and decrease with the sensor's height $z_a$. The tilt modulation modifies the intensity $\sigma$ of the backscattered signal according to the local angle of incidence of the radar beam on the surface. A geometrical approach can also be used to express this intensity (Nieto Borge et al. 2004) and reads

$$ \sigma = c_1 \frac{u \cdot n}{|u| |n|} + c_2, \qquad (2.1) $$

where $n$ and $u$ are the surface-normal and radar-beam vectors (see figure 1), respectively, and $|u|$ ($|n|$) denotes the Euclidean norm of $u$ ($n$). Here, $c_1$ is a scaling factor and $c_2$ an offset, which both depend on the calibration of the radar system and environmental conditions. The determination of these scaling coefficients is an active field of research. To name the most common methods, they can be found based on the signal-to-noise ratio of
<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Typical spatio-temporal distribution of observations from a marine radar scanning irregular unidirectional waves travelling from the left to the right. (b) Random spatio-temporal distribution using the same number of observations.</figcaption></figure>

<table><caption>Table 1. Definition of the three investigated types of observations.</caption><thead><tr><td>Data type</td><td>Observed quantity</td><td>Distribution</td></tr></thead><tbody><tr><td>1</td><td>η</td><td>Random</td></tr><tr><td>2</td><td>η</td><td>Radar like</td></tr><tr><td>3</td><td>σ</td><td>Radar like</td></tr></tbody></table>

the backscattered signal (e.g. Lund et al. 2014), on the shadowing (e.g. Ludeno & Serafino 2019) or statistical (e.g. Gangeskar 2000) properties of radar images or on an external reference measurement (e.g. Naaijen et al. 2018). The determination of these coefficients is outside of the scope of the present paper, hence they are here assumed to be known.

To study the influence of the modulation mechanisms on the wave prediction accuracy, surface reconstructions are systematically performed using three types of wave measurements, referred to as data types. Data type 1 corresponds to a set of randomly distributed surface elevations. Data type 2 refers to surface elevations that are subject to shadowing modulation. Both distributions are depicted in figure 2 with the same number of points. Finally, data type 3 designates radar intensities calculated from (2.1). All data types are summarized in table 1. The influence of the shadowing modulation or of the tilt modulation is studied by comparing results obtained with data types 1 and 2 or 2 and 3, respectively.

### 2.1.2. Linear tilt modulation

The radar intensity model described by (2.1) is related to the surface elevation and slope in a nonlinear manner. Instead of inverting this nonlinear model during the surface reconstruction, a simplified linear model is used. Linearizing the intensity with respect to the surface slope and elevation (see details in Appendix A), (2.1) reduces to

$$ \sigma = c_1 \left( \cos \Theta + \eta_r \sin \Theta - \eta \frac{\sin^2 \Theta}{R} \right) + c_2, \quad (2.2) $$

where the index notation represents the corresponding derivative. Since the surface elevation in (2.2) is multiplied by the factor $\sin^2 \Theta/R \ll \sin \Theta$, the dominant part of $\sigma$ is proportional to the surface radial slope $\eta_r$. Because the slope contains less information than the surface elevation (i.e. loss of an integration constant), using $\sigma$ as observations instead of $\eta$ is expected to impact the performance of the reconstruction methods.

It should be noted that the simplified approach leading to (2.1) and (2.2) does not consider the full complexity of non-coherent radar measurements, whose deterministic
relationship with the surface elevation has yet to be empirically identified. Using real radar data instead of data type 3 is thus expected to have a negative effect on the surface reconstruction accuracy. Nevertheless, algorithms based on this model have been shown to yield sound results using field measurements (Dankert & Rosenthal 2004; Naaijen et al. 2018).

## 2.2. Highly nonlinear surface reconstruction

The spatial and temporal scales of the problem allow us to use the potential flow theory, which assumes that the flow is irrotational and the fluid inviscid and incompressible. The flow is then described by the Laplace equation $\Delta\phi = 0$ in the fluid domain, with $\phi$ the velocity potential, and by a non-penetration condition $\phi_z = 0$ on the seabed, i.e. at $z = -h$ with $h$ the water depth. On the free surface, and using the surface velocity potential $\phi^s(x, y, t) = \phi(x, y, z = \eta, t)$, the dynamic and kinematic free surface boundary conditions (FSBCs) yield (Zakharov 1968)

$$ 
\left. 
\begin{aligned}
\phi_t^s &= -g\eta - \frac{1}{2} |\nabla\phi^s|^2 + \frac{1}{2} (1 + |\nabla\eta|^2) W, \\
\eta_t &= (1 + |\nabla\eta|^2) W - \nabla\phi^s \cdot \nabla\eta,
\end{aligned}
\right\} (2.3)
$$

where $W = \phi_z|_{z=\eta}$ is the vertical velocity on the free surface, $\nabla = \{\partial_x, \partial_y\}^T$ is the horizontal gradient, $t$ is the time and $g$ is the gravitational acceleration.

### 2.2.1. The HOS method

The HOS method is employed for both the generation of the reference ocean surface and its reconstruction. The core of the HOS method is the expression of the potential as a power series of $\eta$ to a given arbitrary order $M \ge 1$, allowing the formulation of $W$ as a Taylor series involving $\eta$ and $\phi^s$ (Dommermuth & Yue 1987; West et al. 1987). Following the order-consistent formulation from West et al. (1987), the HOS FSBCs, expressed on the free surface, read

$$ 
\left. 
\begin{aligned}
\phi_t^s + g\eta &= -\frac{1}{2} |\nabla\phi^s|^2 + \frac{1}{2} \left[ (W^2)_M + |\nabla\eta|^2 (W^2)_{M-2} \right], \\
\eta_t - \phi_z^s &= W_M^{nl} + |\nabla\eta|^2 W_{M-2} - \nabla\phi^s \cdot \nabla\eta,
\end{aligned}
\right\} (2.4)
$$

in which the linear terms have all been moved to the left-hand side, the subscript of the terms involving $W$ denotes the order of expansion and $W_M^{nl} = W_M - W_1 = W_M - \phi_z^s$. The calculations are performed using a pseudo-spectral approach, in which the physical quantities are expressed by means of basis functions that satisfy the Laplace equation, the seabed condition and the horizontal boundary conditions of the computational domain. Here, the latter boundary conditions are assumed periodic so the velocity potential reads

$$ \phi(\mathbf{x}, z, t) = \mathrm{Re} \left\{ \sum_{i=-\infty}^{+\infty} \sum_{j=-\infty}^{+\infty} A_{ij}^{\phi}(t) \exp(i\mathbf{k}_{ij} \cdot \mathbf{x}) \frac{\cosh\left[|\mathbf{k}_{ij}|(z+h)\right]}{\cosh\left(|\mathbf{k}_{ij}|h\right)} \right\}, \quad (2.5) $$

in which $A_{ij}^{\phi}(t)$ are the complex spectral coordinates of $\phi$, $\mathbf{x} = \{x, y\}^T$ and $\mathbf{k}_{ij} = \{k_{xi}, k_{yj}\}^T$, where $k_{xi} = i2\pi/L_x$ and $k_{yj} = j2\pi/L_y$ with $L_x$ and $L_y$ the lengths of the
<figure><img src="image_5.png" /><figcaption>Figure 3. Characteristics of the Cartesian computational domain used to perform the surface reconstruction with HOS-OCM.</figcaption></figure>

computational domain along the x- and y-directions, respectively. The surface elevation has a similar expression without the z-dependent term. In practice, the sums in (2.5) are truncated to finite numbers $N_x$ and $N_y$. This approach allows the spatial-derivative operations to be performed efficiently due to the use of fast Fourier transforms (FFTs).

### 2.2.2. The HOS-observation coupling method

The highly nonlinear surface reconstruction is performed using the HOS-observation coupling method (HOS-OCM), which numerical basis follows that of Desmars et al. (2022). It consists of the inversion of a system that couples the equations describing the evolution of the ocean surface to a set of observation constraints having the form

$$ \mathcal{M}(\phi^s, \eta) = \zeta, \qquad (2.6) $$

in which $\zeta$ is the observed quantity (here, either the surface elevation or the radar intensity), and $\mathcal{M}$ is the observation-mapping function. This function both maps the physical quantities of interest ($\eta$ and $\phi^s$) to the observation space and interpolates them to the spatio-temporal location of the observations. The wave equations are solved on a computational domain that has spatial lengths $L_x$ and $L_y$ and a duration $T_a$ (also called assimilation period). It is discretized with constant steps $\Delta x$ and $\Delta y$ in space and $\Delta t$ in time (see figure 3). It contains $N_x$ and $N_y$ points along the spatial dimensions and $N_t$ points along the temporal dimension, leading to $N_p = N_x N_y N_t$ grid points in total. The reference point in space ($x_0, y_0$) corresponds to the beginning of the computational domain, and the reference point in time $t_0$ corresponds to its first time (see figure 3).

The evolution equations (2.4) are coupled to (2.6) to form the inverse problem that is to be solved for the surface reconstruction. By recasting the resulting system in matrix form, it leads to

$$ \mathbf{A} \cdot \begin{bmatrix} \phi^s \\ \eta \end{bmatrix} = \mathbf{b} \iff \begin{bmatrix} \mathbf{A}^{\partial_t} & g \mathbf{I} \\ -\mathbf{A}^{\partial_z} & \mathbf{A}^{\partial_t} \\ \mathbf{A}^{\mathcal{M}} & \mathbf{I} \end{bmatrix} \cdot \begin{bmatrix} \phi^s \\ \eta \end{bmatrix} = \begin{bmatrix} \mathbf{b}^d \\ \mathbf{b}^k \\ \zeta \end{bmatrix}, \qquad (2.7) $$

where $\mathbf{A}^{\partial_t}$ corresponds to the time-derivative operator matrix, $\mathbf{A}^{\partial_z}$ to the vertical-derivative operator matrix and $\mathbf{I}$ to the identity matrix, all being $N_p$-by-$N_p$ matrices. The $N_o$-by-$2N_p$ (with $N_o$ the number of observations) matrix $\mathbf{A}^{\mathcal{M}}$ corresponds to the observation-mapping
function matrix. Vectors $\phi^s$ and $\eta$ contain the $N_p$ unknown grid values of surface potential and elevation, respectively. The right-hand side terms $b^d$ and $b^k$ correspond to the nonlinear parts of the dynamic and kinematic HOS FSBCs, respectively, and yield

$$b^d = -\frac{1}{2} \left[ (\mathbf{A}^{\partial_x} \cdot \boldsymbol{\phi}^s)^2 + (\mathbf{A}^{\partial_y} \cdot \boldsymbol{\phi}^s)^2 \right] \\ + \frac{1}{2} \left\{ (\mathbf{W}^2)_M + \left[ (\mathbf{A}^{\partial_x} \cdot \boldsymbol{\eta})^2 + (\mathbf{A}^{\partial_y} \cdot \boldsymbol{\eta})^2 \right] (\mathbf{W}^2)_{M-2} \right\}, \quad (2.8)$$

$$b^k = W_M^{nl} + \left[ (\mathbf{A}^{\partial_x} \cdot \boldsymbol{\eta})^2 + (\mathbf{A}^{\partial_y} \cdot \boldsymbol{\eta})^2 \right] W_{M-2} \\ - \left[ (\mathbf{A}^{\partial_x} \cdot \boldsymbol{\eta}) (\mathbf{A}^{\partial_x} \cdot \boldsymbol{\phi}^s) + (\mathbf{A}^{\partial_y} \cdot \boldsymbol{\eta}) (\mathbf{A}^{\partial_y} \cdot \boldsymbol{\phi}^s) \right], \quad (2.9)$$

where $\mathbf{A}^{\partial_x}$ and $\mathbf{A}^{\partial_y}$ are the horizontal-derivative operator matrices, and $W_M$ is a vector whose elements are the vertical velocities $W_M$. Finally, the vector $\boldsymbol{\zeta}$ contains the observations. When not specified by a distinctive mathematical symbol, operations on vectors are performed element-wise. The formulation of the matrices is detailed in Appendix B. The solution vectors $\boldsymbol{\phi}^s$ and $\boldsymbol{\eta}$ are found by minimizing, with the loose generalized minimum residual (LGMRES) algorithm (Baker, Jessup & Manteuffel 2005), the norm of the residual function

$$F = \mathbf{A}^{\mathrm{T}} \cdot \mathbf{A} \cdot \left[ \frac{\boldsymbol{\phi}^s}{\boldsymbol{\eta}} \right] - \mathbf{A}^{\mathrm{T}} \cdot \mathbf{b}. \quad (2.10)$$

For $M=1$, $b^d$ and $b^k$ are full of zeros. For $M > 1$, the solutions are obtained by inverting the system order by order, evaluating $b^d$ and $b^k$ from the lower-order solution. In the following, the first- and third-order solutions, referred to as HOSM1 and HOSM3, respectively, are investigated. The tolerances and restart parameter of the LGMRES iterative process have been fixed in such a way that they lead to converged results in terms of reconstruction accuracy.

Note that the computational cost of HOS-OCM is not related to that of a large dense matrix inversion as one might initially assume. Indeed, the chosen iterative solver for the minimization of (2.10) allows matrix-free calculations: **A** is not assembled, only matrix-vector products involving its submatrices are performed. These submatrices exhibit a high level of sparsity that leads them to have a product cost $O(N_p N_s)$ for derivation or $O(N_o N_s)$ for interpolation with $N_s \ll N_p$ the number of stencils to perform the considered operation. Matrix-vector products involving the spatial-derivative operator matrices $\mathbf{A}^{\partial_{x,y,z}}$ are accelerated with FFTs resulting in a cost $O(N_p \log N_s)$. The transpose operations require the exact same computational effort. In fact, there is no additional cost coming from the matrix formulation of the problem, and the resulting number of operations is of the same order of magnitude as that of other HOS-based prediction algorithms relying on variational assimilation.

However, in contrast to variational approaches that control the HOS variables (namely the surface elevation and potential) at only one time step and propagate the solution over the assimilation period with a classical time-stepping procedure (typically using an explicit fourth-order Runge-Kutta scheme), HOS-OCM controls the variables at all time steps simultaneously. This formulation has two interesting numerical properties. First, no sequential time propagation of the solution is required, making possible the time parallelization of the model evaluation in the optimization process. This parallelization in
time offers a greater impact on the computational time reduction than the parallelization in space whose effect is limited by the relatively small domain used for such surface reconstructions. Second, the embedded time integration of the wave equations has the properties of implicit schemes at no additional cost. More specifically, it exhibits high numerical stability – no evidence of instability was observed even for very large time steps – without the need to perform the additional system inversion that is required when the solution is sequentially propagated with implicit schemes. This gives flexibility in the choice of the time-step size, thus in the management of the computational time.

## 2.3. Parameterization of analytical wave models

The other investigated surface reconstruction method relies on the parameterization of analytical wave models. Besides a linear model that represents the water surface as the superposition of independent sine waves subject to the linear dispersion relation, two weakly nonlinear models are investigated. Based on a Lagrangian description of the surface motion, the improved choppy wave model catches nonlinear effects on both shape, i.e. sharper crests and flatter troughs, and phase, i.e. higher wave velocity (Guérin et al. 2019). In contrast, the second weakly nonlinear model considered here includes the correction of the dispersion relation pertaining to the improved choppy wave model, but no shape correction.

For these models, detailed below, the surface reconstruction is based on a least-squares algorithm that minimizes a quadratic cost function, similar to the approaches investigated by Desmars et al. (2020). The cost function that is minimized reads

$$C = \frac{1}{2} |\mathbf{\Lambda}(\mathbf{p}) - \mathbf{\Lambda}^o|^2, \qquad (2.11)$$

where the vector $\mathbf{\Lambda}$ contains the wave model estimates of the observations, $\mathbf{p}$ contains the parameters of the wave model and $\mathbf{\Lambda}^o$ contains either observed elevations ($\eta^o$) or radar intensities ($\sigma^o$) depending on the chosen data type. The optimal wave model parameters are obtained by applying the Levenberg-Marquardt algorithm (Moré 1978) with explicit formulations of the Jacobian matrices that allow the procedure to be performed efficiently. The Jacobian matrices depend on the considered wave model and data type. Their expression for each configuration is given in Appendix C. In the rest of the paper, and in contrast to HOS-OCM, this reconstruction algorithm is referred to as the analytical-model-based method (AMBM).

### 2.3.1. Linear wave theory

The linear wave theory (LWT) gives the surface velocity potential and corresponding elevation according to

$$\phi_{LWT}^{s}(\mathbf{x}, t) = \sum_{n=1}^{N} \frac{\omega_n}{k_n} (a_n \sin \psi_n - b_n \cos \psi_n), \qquad (2.12)$$

$$\eta_{LWT}(\mathbf{x}, t) = \sum_{n=1}^{N} (a_n \cos \psi_n + b_n \sin \psi_n), \qquad (2.13)$$

in which $\psi_n = \mathbf{k}_n \cdot \mathbf{x} - \omega_n t$, with $\mathbf{k}_n = \{k_{xn}, k_{yn}\}^T$, and $a_n$ and $b_n$ are the model parameters related to the amplitude and phase of N individual wave components of predefined wavenumbers $k_n = |\mathbf{k}_n|$ related to angular frequencies $\omega_n$ through the linear dispersion relation $\omega_n^2 = gk_n \tanh(k_n h)$.
### 2.3.2. Improved choppy wave model

The improved choppy wave model (ICWM) was developed and characterized by Guérin et al. (2019). It relies on the Lagrangian description of the water surface which allows the derivation, in a relatively simple mathematical form, of terms appearing at higher orders of expansion in wave steepness in the Eulerian formalism. Compared with the second-order Lagrangian solution (e.g. Nouguier, Chapron & Guérin 2015), ICWM corrects the nonlinear effects affecting the waves’ velocity resulting in a more accurate description of the free surface kinematics. An Eulerian form of ICWM was used by Desmars et al. (2020) to predict unidirectional waves and writes

$$ \phi_{\text{ICWM}}^{s}(\mathbf{x}, t) = \sum_{n=1}^{N} \frac{\tilde{\omega}_{n}^{s}}{k_{n}} (a_{n} \sin \Psi_{n} - b_{n} \cos \Psi_{n}), \quad (2.14) $$

$$ \eta_{\text{ICWM}}(\mathbf{x}, t) = \sum_{n=1}^{N} (a_n \cos \Psi_n + b_n \sin \Psi_n) + \frac{1}{2} \sum_{n=1}^{N} (a_n^2 + b_n^2) k_n, \quad (2.15) $$

with

$$ \left. \begin{aligned} \Psi_n &= \mathbf{k}_n \cdot \left[ \mathbf{x} - \sum_{m=1}^{N} \frac{\mathbf{k}_m}{k_m} \left( -a_m \sin \psi_m^{nl} + b_m \cos \psi_m^{nl} \right) \right] - \omega_n^{nl} t, \\ \psi_n^{nl} &= \mathbf{k}_n \cdot \mathbf{x} - \omega_n^{nl} t, \\ \omega_n^{nl} &= \omega_n + \frac{1}{2} \mathbf{k}_n \cdot \sum_{m=1}^{N} (a_m^2 + b_m^2) \omega_m \mathbf{k}_m, \\ \tilde{\omega}_n^s &= \omega_n - \frac{1}{2} \mathbf{k}_n \cdot \sum_{m=1}^{N} (a_m^2 + b_m^2) \omega_m \mathbf{k}_m. \end{aligned} \right\} \quad (2.16) $$

Note that the formulation of the surface potential of ICWM, not given yet in the literature, is discussed in Appendix D.

### 2.3.3. Linear wave theory with corrected dispersion relation

To quantify the importance of the correction of the dispersion relation only, a linear wave model that includes the correction of the dispersion relation pertaining to ICWM is also used. This model is referred to as LWT with corrected dispersion relation (LWT-CDR) and reads

$$ \phi_{\text{LC}}^{s}(\mathbf{x}, t) = \sum_{n=1}^{N} \frac{\tilde{\omega}_{n}^{s}}{k_{n}} (a_{n} \sin \psi_{n}^{nl} - b_{n} \cos \psi_{n}^{nl}), \quad (2.17) $$

$$ \eta_{\text{LC}}(\mathbf{x}, t) = \sum_{n=1}^{N} (a_n \cos \psi_n^{nl} + b_n \sin \psi_n^{nl}). \quad (2.18) $$

## 2.4. Numerical differences between HOS-OCM and AMBM

Throughout the paper, a systematic comparison of the reconstruction and prediction accuracy of each algorithm is performed. Leaving aside hydrodynamic considerations, this
## N. Desmars and others

section aims at highlighting the numerical differences between HOS-OCM and AMBM to support the discussion of the results.

Even if both HOS-OCM and AMBM use least squares, through the minimization of $|F|$ (2.10) for HOS-OCM and of $C$ (2.11) for AMBM, a fundamental difference between the two methods is the nature of the control parameters of the minimization process. While they consist of the surface elevation and potential at the grid points for HOS-OCM, AMBM solves for the amplitude and phase of the wave components, and the surface fields are calculated from the model expressions after the minimization. This way, AMBM always gives a solution to the equations of the surface dynamics according to the considered wave model (i.e. LWT, LWT-CDR or ICWM). In contrast, the wave equations for HOS-OCM are seen as constraints that are not necessarily exactly verified by the solution. In addition, the wavenumbers of the wave components are not defined the same way for each method: they are implicitly defined by the grid size and resolution for HOS-OCM, and arbitrarily chosen for AMBM. To improve the comparability of the methods, the wave components of AMBM are chosen according to those implicitly defined by the computational grid of HOS-OCM. Although this might not correspond to the optimal choice of parameters for AMBM, convergence studies on the number of wave components have shown that, for both the synthetic and experimental configurations depicted in this paper (§§ 3 and 4, respectively), this choice leads to results that fairly represent the accuracy of the AMBM solver. Because numerical instabilities are arising for AMBM when very long modes are solved for, wave components whose wavenumber $k$ is lower than half the peak wavenumber $k_p$ of the observed wave field (i.e. the wavenumber of maximal energy) are not considered in the LWT, LWT-CDR and ICWM solutions.

### 2.5. Influence of the spatial periodicity assumption

Both HOS-OCM and AMBM (with similar wavenumbers to those of HOS-OCM, see the previous section) assume that the solution of the surface reconstruction problem is $L_x$- and $L_y$-periodic, while the observations have no reason to exhibit spatial periodicity. Forcing the solution to match at the spatial boundaries of the reconstructed domain leads to discrepancies if no numerical treatment is employed. In the presented applications, the discrepancies induced by the periodicity assumption are restricted to a relatively small region on the domain boundaries. In consequence, no special treatment for the periodicity assumption was used *a priori*, the influence of these discrepancies on the results is reduced by evaluating the quality of the solution over a restricted domain that ends one peak wavelength ($\lambda_p = 2\pi/k_p$) away from each boundary. The size of the regions to be removed from the analysis was chosen for this study based on the evaluation of the actual results of the investigated cases. Changing the configuration, in particular the space and time extent of the observed domain, could modify the size of the impacted regions by the periodicity assumption.

## 3. Reconstructions

In this section, the numerical set-up to evaluate the accuracy of the reconstruction methods from synthetic surface observations is described before discussing the results.

### 3.1. *Numerical set-up*

Two unidirectional wave fields, made of regular and irregular waves that propagate along the positive x-direction, aim at being reconstructed. Throughout the paper, the
reconstruction and prediction configurations are set up with respect to the characteristics of the irregular wave field, even in the case of regular waves that is seen as a ‘simplified’ wave field of similar scales to help interpret the results. The irregular wave field is based on a JONSWAP energy distribution (Hasselmann et al. 1973) with a peak period $T_p = 10$ s and a peak-enhancement factor of 3.3. The wave height varies depending on the investigated wave steepness $H_s/\lambda_p$ (differing by a factor of $\pi$ from the alternative steepness definition $\frac{1}{2}k_pH_s$), with $H_s$ the significant wave height defined as four times the standard deviation of the surface elevation. The regular wave field is defined from the peak wave characteristics of the irregular one, i.e. the period is $T = T_p$ and different values of wave steepness $H/\lambda (=kA/\pi)$ are investigated, with $H = 2A$ the wave height defined as the crest-to-trough height, and $\lambda = 2\pi/k = \lambda_p$ the wavelength. These wave fields, called reference wave fields and indicated by the superscript ‘r’ (e.g. $\eta^r$), are generated using the HOS solver HOS-ocean (Ducrozet et al. 2016) with a nonlinear order $M=4$, which provides converged results in terms of nonlinearity. The length of the domain is set to 32 peak wavelengths, deep water is assumed and a relaxation period $T_r$, such as described by Dommermuth (2000), is used to smoothly turn the linearly initialized wave field into a fourth-order solution. A conservative simulation time of $2T_r$ is used before considering that the nonlinearities are fully developed and generating the synthetic observations. Regarding the spatial discretization, 32 points per peak wavelength are used, which ensures that the shortest wave sampled by the radar is properly resolved. This discretization, however, limits the steepness of the simulated irregular wave field to $H_s/\lambda_p \lesssim 3.2\%$ according to application ranges of highly nonlinear potential flow solvers (Ducrozet, Bonnefoy & Perignon 2017).

To generate the radar observations, a virtual radar scans the surface and gives at every non-shadowed point one of the two observed quantities, that is either the surface elevation $\eta^o = \eta_j^r$ (with $j = 1, \dots, N_o$) or the intensity $\sigma^o = \sigma_j^r$ according to (2.1). In the case the observations consist of radar intensities, the distances $R^o = R_j$ and angles of incidence $\Theta^o = \Theta_j$ (see figure 1) are stored as well. In operating conditions, $R$ is approximated based on the time of flight of the radar beam, and $\Theta$ is calculated based on the geometric relation $\Theta = \cos^{-1}(z_a/R)$. Moreover, to decouple the effects of the wave shadowing (i.e. larger gaps in the spatial distribution of observations) from those of the wave steepness (i.e. stronger impact of the nonlinear wave physics) on the reconstruction accuracy for larger wave heights, the height of the radar antenna is made proportional to the wave height according to $z_a = 8H_s$ in the case of (ir)regular waves. This way, the shadowing properties are kept similar for all the investigated wave heights, and the observed differences in reconstruction accuracy between various steepness values are only related to hydrodynamic nonlinearities. The influence of the wave shadowing is quantified by comparing the results with those obtained with randomly distributed observations (i.e. data type 1 vs data type 2). Since the number of observations in the case of data type 1 is set similar to that obtained with the radar sampling of the surface (which includes wave shadowing), the overall density of observations is independent of the data type. Despite the unrealistic nature of the dependence between $z_a$ and the wave height, the factor of 8 already depicts a situation of strong shadowing (see figure 2a). The radar horizontal location is chosen such that it is facing the incoming waves and that there is a gap of $3\lambda_p \approx 470$ m between the radar and the reconstructed domain, i.e. $x_a = x_0 + L_x + 3\lambda_p$. The remaining radar characteristics are chosen according to Naaijen & Wijaya (2014), i.e. the range resolution is set to 7.5 m and the antenna rotation period, fixing the temporal resolution, to 1.5 s. Only the observations that fall into the reconstructed domain $[x_0, x_0 + L_x] \times [t_0, t_0 + T_a]$ are used for the reconstruction.
## N. Desmars and others

Concerning the computational domain of HOS-OCM, its spatial extent is $L_x = 12\lambda_p$ for a number of points $N_x = 246$, leading to a resolution of $\Delta x \approx \lambda_p/20.4 \approx 7.6$ m. This is in agreement with the maximal spatial resolution of 7.5 m imposed by the observations in the case the ocean surface is sampled by the radar. The assimilation period is $T_a = 3T_p$, which is long enough to give a converged reconstruction accuracy for the chosen measurement length, and the number of time steps is $N_t = 48T_a/T_p = 144$ (i.e. $\Delta t \approx 0.21$ s).

## 3.2. Reconstruction accuracy quantification

For the quantification of the error between the reconstructed surface and the reference, the surface similarity parameter (SSP) is used. This indicator was developed by Perlin & Bustamante (2016) and already used in the context of ocean wave prediction (e.g. Lünser et al. 2022). For the surface elevation, the SSP between the reconstructed field $\eta$ and the reference $\eta^r$ is formulated as

$$ \text{SSP}_{\eta}(t) = \frac{\left( \int |\mathcal{F}_{\eta}(k, t) - \mathcal{F}_{\eta^r}(k, t)|^2 dk \right)^{1/2}}{\left( \int |\mathcal{F}_{\eta}(k, t)|^2 dk \right)^{1/2} + \left( \int |\mathcal{F}_{\eta^r}(k, t)|^2 dk \right)^{1/2}}, \quad (3.1) $$

where $\mathcal{F}_{\eta}$ denotes the spatial FFT of $\eta(x, t)$ over the points of the computational domain of HOS-OCM. A similar expression holds for the surface potential error SSP$_{\phi^s}(t)$. As mentioned in § 2.5, a restricted region that begins (ends) one peak wavelength after (before) the left (right) spatial boundary is selected for the quantification of the reconstruction accuracy to limit the influence of the periodicity assumption, i.e. $x \in [x_0 + \lambda_p, x_0 + L_x - \lambda_p]$. The value of the SSP is bounded by 0 and 1, meaning a perfect agreement or disagreement between the compared signals, respectively. Because the motion response of marine structures is most likely influenced by a limited frequency range, i.e. similar to a low-pass filter, the presented SSP calculation focuses on wavenumbers whose expected impact is significant. More specifically, this is done by retaining in (3.1) only the complex amplitudes $\mathcal{F}$ associated with wavenumbers lower than or equal to a cutoff value, here defined as $5k_p$. This choice is supported by the fact that the wave field evolution is mainly driven by wave components of wavenumber $k < 5k_p$ (Ducrozet et al. 2017). Moreover, to have an error indicator representing the global accuracy over the reconstructed domain, the time-averaged SSP between the two instants $t_0$ and $t_0 + T_a$ is calculated. For the surface elevation, it reads

$$ \overline{\text{SSP}}_{\eta} = \frac{1}{T_a} \int_{t_0}^{t_0+T_a} \text{SSP}_{\eta}(t) \, dt. \quad (3.2) $$

An expression similar to (3.2) is used to quantify $\overline{\text{SSP}}_{\phi^s}$. In the case of irregular waves, the SSP values are averaged over 200 wave field realizations with different sets of randomly chosen initial wave phases.

## 3.3. Reconstruction results

In this section, the ability of HOS-OCM and AMBM to accurately reconstruct the observed wave field is studied. First, regular waves are considered, leading to a clear distinction between the characteristics of each modelling method. Then irregular waves are reconstructed, showing the performance of the methods in more realistic wave conditions.
### 3.3.1. Regular waves

Reconstruction results for regular waves are presented in figure 4. They show that, for all data types, the reconstruction accuracy decreases for larger wave steepness when the wave model is not able to catch enough nonlinear features, consisting here of increasing the wave shape asymmetry and the wave velocity, the latter appearing from the third order. The HOSM1 and LWT are fully linear and miss the aforementioned nonlinear effects, penalizing the reconstruction accuracy when the steepness increases. The trend followed by LWT-CDR is almost similar to that of LWT (and of HOSM1), showing that modelling the nonlinear correction of the wave velocity is not enough to improve the reconstruction significantly. However, the influence of the wave velocity having linear and quadratic dependencies on the simulation time and wave steepness, respectively, means that the nonlinear correction of LWT-CDR is expected to be more impactful for longer assimilation period $T_a$ and larger steepness. In contrast, the accuracy is very high for HOSM3 and ICWM, and barely affected by the steepness variation. This is because both models include all the nonlinear effects associated with regular waves up to the third order. As seen from the slightly lower accuracy obtained with data type 2 compared with data type 1, the wave shadowing effect, which tends to concentrate the discrepancies at the non-visible parts of the waves (e.g. troughs), has a limited impact. Using observations of radar intensity (data type 3) instead of elevation (data type 2) has a small influence as well, which validates that both HOS-OCM and AMBM can extract correct information about the surface dynamics from elevation or slope samples. In the case of data type 3, the fact that the accuracy of HOS-OCM is generally a little lower than that of AMBM most likely comes from the numerical differences detailed in § 2.4 between the two methods. In particular, performing the reconstruction from radar intensities instead of elevations might affect the relative weight of the HOS equations (2.4) with respect to the observation-mapping function in the minimized residual (2.10) of HOS-OCM, resulting in a lower accuracy. A penalty approach (e.g. Nocedal & Wright 1999) could be employed to investigate the effect of the wave model enforcement (i.e. forcing the solution to strictly satisfy the HOS equations) on the reconstructed surface. Nevertheless, the SSP values remain below 0.05 in all configurations, and almost the same results are obtained for the surface potential.

Reconstruction of regular waves of larger steepness (up to $H/\lambda \sim 10\%$) have been performed and provide results that are in line with the physical properties of the wave models, i.e. similar accuracy between HOSM3 and ICWM and progressive distinction between LWT and LWT-CDR for increasing steepness. For the very steep cases, however, the reference data generation, surface reconstruction and surface prediction procedures had to be slightly adapted to prevent various numerical perturbations from appearing, and whose study is left out of this paper for the sake of brevity.

### 3.3.2. Irregular waves

Reconstruction results for irregular waves are presented in figure 5. As expected, and similar to the case of regular waves, the lower-order models are less accurate, and the SSP values increase for larger wave steepness. However, larger SSP values than those of the regular waves are found, which can be explained by the increased complexity resulting from the interaction phenomena between waves of different scales. This is especially clear when the reconstruction relies on slope information (i.e. data type 3) for which the accuracy strongly depends on the correct resolution of the modulation of short waves by long waves, making the whole range of wave scales important to resolve. In contrast, mainly the waves of high amplitude (localized around the peak of the energy spectrum)
<figure><img src="image_6.png" /><figcaption>Figure 4. Reconstruction accuracy for the regular wave field as a function of the steepness $H/\lambda$ for HOSM1 (---●--), HOSM3 (—●—), LWT (-Δ-·, red), LWT-CDR (··Δ··, red) and ICWM (-·Δ-·, red). Results using data type (a,d) 1, (b,e) 2 and (c,f) 3 are presented (see table 1 for definitions).</figcaption></figure>

influence the correct extraction of wave information from the observations when they consist of surface elevations (i.e. data types 1 and 2), simplifying the reconstruction. Radar intensities are also harder to use for surface reconstruction because they contain less information about the wave field than surface elevations, as explained in § 2.1.2. Regarding the linear approaches, HOSM1 systematically leads to a higher accuracy than LWT for data types 1 and 2, and to a significantly lower accuracy for data type 3. This could be explained, as for the regular waves, by the numerical differences between HOS-OCM and AMBM (see § 2.4). In combination with the numerical characteristics of HOS-OCM, the third-order modelling of HOSM3 gives the most accurate reconstructions for data types 1 and 2. As shown in figure 6(a) for data type 2, the improvement of HOSM3 over HOSM1 is mainly visible for large surface deformation. Without providing clear hints for interpretation, the discrepancies implied by the use of data type 3 in nonlinear reconstructions of the same surface can be spotted in figure 6(b). In general, the SSP values for the surface potential are a bit larger than for the surface elevation, which could come from the fact that, for all data types, the observations are related to the surface elevation but not directly to the surface potential. The potential is calculated from the wave model only, leading its reconstruction error to be slightly higher than that of the elevation.

## 4. Predictions

After investigating the performance of the methods for surface reconstruction, this section focuses on quantifying the accuracy that can be obtained by propagating the reconstructed solutions in view of predictions.

### 4.1. *Definition of the prediction domain*

The prediction accuracy is quantified over the theoretically accessible prediction zone corresponding to the spatio-temporal domain within which pieces of information about the
Nonlinear prediction of remotely measured ocean waves

<figure><img src="image_7.png" /><figcaption>Figure 5. Reconstruction accuracy for the irregular wave field as a function of the steepness $H_s/\lambda_p$ for HOSM1 (---), HOSM3 (●-●-●), LWT (---Δ--), red), LWT-CDR (···Δ···, red) and ICWM (·-Δ-·, red). Results using data type (a,d) 1, (b,e) 2 and (c,f) 3 are presented (see table 1 for definitions). The shaded areas show the standard deviation for HOSM3 and ICWM for the 200 wave field realizations.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Reconstructed surface elevation of irregular waves of steepness $H_s/\lambda_p \approx 3.2\%$ at time $t_0 + T_a$ using observations subject to wave shadowing and that consist of (a) surface elevations (i.e. data type 2) and (b) radar intensities (i.e. data type 3). Comparison of the reference surface (—, green) with HOSM1 (---, SSP$_\eta \approx 0.078$), HOSM3 (—, SSP$_\eta \approx 0.024$ (a) and 0.076 (b)) and ICWM (—, red, SSP$_\eta \approx 0.062$).</figcaption></figure>

reconstructed wave components, travelling at the group velocity $c_g = \partial\omega/\partial k$, overlap (e.g. Naaijen, Trulsen & Blondel-Couprie 2014; Qi et al. 2018b). In practice, the boundaries of the prediction zone are defined by considering that the wave field can be accurately described by a finite frequency bandwidth containing the most energetic wave components. Here, similar to Desmars et al. (2020) and Huchet et al. (2021), cutoff angular frequencies $\omega_\ell$ and $\omega_h$ (with $\omega_\ell < \omega_h$) are selected such that $E(\omega_\ell) = E(\omega_h) = 0.05E(\omega_p)$, with $E(\omega)$ the wave energy spectrum of the irregular wave field and $\omega_p = 2\pi/T_p$. These angular frequencies give the highest and lowest group velocities of the wave components of interest that define the prediction zone boundaries. Because it is supposed that the
<figure><img src="image_9.png" /><figcaption>Figure 7. Definition of the theoretically accessible prediction zone from information contained in the reconstructed solution at $t = t_0 + T_a$.</figcaption></figure>

radar is mounted close to the location where the waves have to be anticipated, the target location for prediction is the radar location $x_a$ and the predicted surface is evaluated as a time series. Also, to maximize the prediction horizon, the last reconstructed solution (i.e. at $t = t_0 + T_a$) is taken as the initial condition for the propagation. Finally, since the solution is considered valid over a restricted domain due to the spatial periodicity assumption (see § 2.5), the prediction zone is determined by neglecting information out of $[x_0 - \lambda_p, x_0 + L_x + \lambda_p]$. In consequence, a future time $t$ is included in the prediction zone at $x_a$ if it is part of the interval $[t_h, t_\ell]$ with

$$t_h = t_0 + T_a + \frac{x_a - (x_0 + L_x - \lambda_p)}{c_{gh}}, \quad (4.1)$$

$$t_\ell = t_0 + T_a + \frac{x_a - (x_0 + \lambda_p)}{c_{gl}}, \quad (4.2)$$

and the group velocities are calculated using the linear dispersion relation, i.e. $c_g = \frac{1}{2}\omega/k[1 + 2kh/\sinh(2kh)]$, which already gives a good estimation of the extent of the prediction zone compared with a higher-order formulation of the dispersion relation (Qi et al. 2018a). Figure 7 shows the geometrical interpretation of the prediction zone definition. In the depicted configuration, the obtained limits of the prediction horizon lead to $t_h - (t_0 + T_a) \approx 14T_p$ and $t_\ell - t_h \approx 6T_p$.

## 4.2. Numerical aspects of the prediction using HOS-OCM

While the analytical models can be propagated in space and time without requiring any special numerical treatment, the HOS solutions from HOS-OCM need a fixed, periodic spatial domain to be propagated over. Here, the propagation of the HOS solutions is done as follows. First, the spatial computational domain of HOS-OCM is extended to include the location of the prediction $x_a$. This domain is set similarly to that used for the generation of the reference surface described in § 3.1, and it is used to express the reconstructed solution at $t = t_0 + T_a$ with zeros for $x \notin [x_0, x_0 + L_x]$. The prediction is thus based on the propagation of the reconstructed solution obtained at the end of the assimilation period. Second, to avoid discontinuities, the solution is made smoothly go to zero at both ends of
the reconstructed domain. More specifically, the surface fields $\eta$ and $\phi^s$ are multiplied by the third-order polynomial

$$f_s(x) = \begin{cases} \left[ \frac{x - (x_0 - \lambda_p)}{\lambda_p} \right]^2 \left\{ 3 - 2 \left[ \frac{x - (x_0 - \lambda_p)}{\lambda_p} \right] \right\}, & \text{if } x \in [x_0 - \lambda_p, x_0], \\ 1, & \text{if } x \in [x_0, x_0 + L_x], \\ 1 - \left[ \frac{x - (x_0 + L_x)}{\lambda_p} \right]^2 \left\{ 3 - 2 \left[ \frac{x - (x_0 + L_x)}{\lambda_p} \right] \right\}, & \text{if } x \in [x_0 + L_x, x_0 + L_x + \lambda_p]. \end{cases} \quad (4.3)$$

Finally, the HOS-ocean solver is used for the forward propagation of the solution with an order of nonlinearity $M$ matching that used for the reconstruction, i.e. either 1 or 3.

As described in § 4.1, even if the initial surface is perfectly known over a limited spatial domain, the dispersion of wave information during the surface propagation automatically prevents the prediction at a further point in space from perfectly matching the true solution. The accuracy of such a prediction depends on the amount of truncated energy by the cutoff frequencies defining the extent of the prediction zone, and a perfect prediction remains out of reach even with a perfectly reconstructed surface. Consequently, the optimal prediction accuracy can be quantified by applying the propagation procedure described above (although with $M=4$) to the reference solution and comparing the obtained prediction with the normally propagated reference solution. The determination of this error is done in the following for the case of irregular waves and indicates how far the accuracy obtained with the investigated prediction algorithms is from the best achievable one. This optimal prediction accuracy quantification does not apply to regular wave fields because they are not subject to wave dispersion.

Because propagating a solution which boundaries have been smoothed down to zero potentially introduces non-physical artefacts in the simulation, the spectra of the reference prediction and that of the optimal prediction have been compared. This comparison showed that no artefact appears in the solution subject to the described smoothing process.

### 4.3. Prediction accuracy quantification

The SSP is used again to quantify the accuracy of the predicted surface elevation and potential. Because the comparison with the reference surface relies on time series, its formulation is different from (3.1) and follows

$$\widetilde{\text{SSP}}_{\eta} = \frac{\left( \int |\mathcal{F}_{\eta}(f) - \mathcal{F}_{\eta^r}(f)|^2 \, df \right)^{1/2}}{\left( \int |\mathcal{F}_{\eta}(f)|^2 \, df \right)^{1/2} + \left( \int |\mathcal{F}_{\eta^r}(f)|^2 \, df \right)^{1/2}}, \quad (4.4)$$

where $\mathcal{F}_{\eta}$ denotes the temporal FFT of $\eta(x_a, t)$ over the prediction period, i.e. $t \in [t_h, t_l]$, and $f$ is the frequency. A similar expression holds for the surface potential error $\widetilde{\text{SSP}}_{\phi^s}$. To only take into account the theoretically predictable wave components, the integral boundaries span the frequency range that is used to determine the extent of the prediction zone, i.e. the interval between $f_l = \omega_l / (2\pi)$ and $f_h = \omega_h / (2\pi)$.
N. Desmars and others

<figure><img src="image_10.png" /><figcaption>Figure 8. Prediction accuracy for the regular wave field as a function of the steepness $H/\lambda$ for HOSM1 (---), HOSM3 (—●—), LWT (----Δ—), red), LWT-CDR (···Δ···, red) and ICWM (··Δ·, red). Results using data type (a,d) 1, (b,e) 2 and (c,f) 3 are presented (see table 1 for definitions).</figcaption></figure>

4.4. Prediction results

Similar to the reconstruction, the ability of the different methods to predict the observed
wave field is investigated for both regular and irregular wave fields.

4.4.1. *Regular waves*

Prediction results for regular waves are presented in figure 8, and unambiguously show the
advantage of models that take into account the nonlinear correction of the wave velocity.
Indeed, the accuracy decreases significantly with the wave steepness for both HOSM1
and LWT, while it remains very high for all the other models. This is in agreement with
previous studies that compare the prediction performance of linear and nonlinear models
and find that correct predictions of severe sea states can only be obtained with an accurate
modelling of the wave velocity (e.g. Blondel-Couprie et al. 2013; Desmars et al. 2020;
Meisner et al. 2023). The effect of the nonlinear phase shift on the predicted surfaces is
depicted in figure 9. The influence of both measurement modulation mechanisms, i.e. the
wave shadowing effect and the use of radar intensities, on the prediction accuracy is very
limited.

4.4.2. *Irregular waves*

Prediction results for irregular waves are presented in figure 10. In this configuration,
the benefit of HOSM3 is clear. Even if LWT-CDR and ICWM significantly improve the
prediction compared with HOSM1 and LWT, the additional modelling effects included
in HOSM3 substantially increase the accuracy. For data types 1 and 2, the results
obtained with HOSM3 are satisfactory (see surface profiles in figure 11a), although
there remains room to reach the optimal accuracy defined in § 4.2. Even for data
type 3, for which HOSM3 leads to a relatively poor reconstruction compared with
ICWM (see figure 5c,f), the propagation performance of HOSM3 leads to a prediction
accuracy comparable to or even higher than that of ICWM for large wave steepness.
Nonlinear prediction of remotely measured ocean waves

<figure><img src="image_11.png" /><figcaption>Figure 9. Predicted surface elevation of regular waves of steepness $H/\lambda \approx 3.2\%$ at location $x_a$ from a reconstruction using observations subject to wave shadowing and that consist of radar intensities (i.e. data type 3). (a) Comparison of the reference surface (—, green) with LWT (---, red, $\overline{SSP}_\eta \approx 0.287$) and LWT-CDR (······, red, $\overline{SSP}_\eta \approx 0.049$). (b) Comparison of the reference surface with HOSM3 (—, $\overline{SSP}_\eta \approx 0.014$) and ICWM (----, red, $\overline{SSP}_\eta \approx 0.023$).</figcaption></figure>

Examples of nonlinear predictions using data type 3 are shown in figure 11(b) to illustrate
the impact of the complete third-order modelling compared with ICWM. The nonlinear
phase velocity correction of ICWM is independent of the properties of individual wave
components, i.e. it acts as a global horizontal shift of the surface according to the sum of
every wave component’s Stokes drift (Guérin *et al.* 2019). In consequence, ICWM only
models an ‘average’ nonlinear phase shift that is mainly driven by the most energetic
wave components, and misses the frequency-dependent phase correction modelled by
HOSM3. More specifically, compared with the modification of the dispersion relation
due to third-order resonant effects (Longuet-Higgins & Phillips 1962), ICWM tends
to overestimate (underestimate) the velocity correction of low-(high-)frequency wave
components. It is also interesting to note that in the case of data type 3, ICWM gives
better prediction results than LWT-CDR, while they are similar for data types 1 and 2.
This emphasizes the need for accurate modelling of surface slopes when the observations
are radar intensities.

**5. Application to experimental data**

In this section, the surface reconstruction/prediction algorithms are applied to real wave
measurements from wave tank experiments. A description and characterization of the
experimental data are provided by Clauss & Klein (2011) and briefly recalled below.

5.1. *Experimental set-up*

The experimental campaign was carried out in the seakeeping wave tank of the Technical
University of Berlin and aimed at studying the spatial development of an extreme wave. To
this end, the New Year wave (also known as the Draupner wave) was reproduced at a scale
1:70 in the wave tank of length 120 m, width 8 m and depth 1 m. By generating the same
unidirectional wave field multiple times and successively moving an array of wave gauges,
a total number of 520 gauge measurements were performed with a spatial step of 0.1 m.
In full scale (which will be used hereafter), the spatial extent of the measurements
<figure><img src="image_12.png" /><figcaption>Figure 10. Prediction accuracy for the irregular wave field as a function of the steepness $H_s/\lambda_p$ for HOSM1 (---), HOSM3 (●-●-●), LWT (---△-·, red), LWT-CDR (···△··, red) and ICWM (·-△-·, red). The optimal accuracy defined in § 4.2 is also shown (—). Results using data type (a,d) 1, (b,e) 2 and (c,f) 3 are presented (see table 1 for definitions). The shaded areas show the standard deviation for HOSM3 and ICWM for the 200 wave field realizations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Predicted surface elevation of irregular waves of steepness $H_s/\lambda_p \approx 3.2\%$ at location $x_a$ from a reconstruction using observations subject to wave shadowing and that consist of (a) surface elevations (i.e. data type 2) and (b) radar intensities (i.e. data type 3). Comparison of the reference surface (—, green) with HOSM1 (---, $\overline{SSP}_\eta \approx 0.498$), HOSM3 (—, $\overline{SSP}_\eta \approx 0.036$ (a) and 0.062 (b)) and ICWM (—, $\overline{SSP}_\eta \approx 0.193$). The optimal prediction defined in § 4.2 is also shown in (a) (—, $\overline{SSP}_\eta \approx 0.010$).</figcaption></figure>

is 3633 m with a resolution of $\Delta x_e = 7$ m, which is comparable to the spatial scale and resolution of radar measurements. The New Year wave appears at predefined target location and time, and the energy distribution of the surrounding wave field is assumed to be described by a JONSWAP spectrum with $T_p = 16.7$ s, $H_s = 11.92$ m and a peak-enhancement factor of 1. With a water depth of 70 m, it leads to a peak wavelength $\lambda_p \approx 364$ m.
Nonlinear prediction of remotely measured ocean waves

This experimental configuration generalizes the applicability of the reconstruction/
prediction algorithms in three ways. First, the dimensionless wavenumber $k_p h \approx 1.2$
makes the configuration fall into the finite water depth regime. Second, even if no
specific characterization of the signal perturbations was done, the data generation process,
based on the (generally non-perfect) repeatability of model tests (e.g. calibration and
measurement inaccuracies with respect to vertical surface displacement but also horizontal
position of the wave gauge array and distance in between) contaminates the space series of
surface elevation with spurious fluctuations. Because the surface elevation was measured
with an array of 13 wave gauges separated by an interval of $2\Delta x_e$, that was translated
by $\Delta x_e$ between two successive wave generations, the main measurement repeatability
errors take the form of localized fluctuations of wavelength $2\Delta x_e$. The amplitude of these
fluctuations depends on the repeatability between the two successive measurements for
which the array of wave gauges was translated by $\Delta x_e$, it is thus not constant all over
the covered spatial domain. For instance, such fluctuations are clearly identifiable in
figure 12 around $x - x_0 = 2800$ m. To check that the presented approaches are applicable
to non-perfect input data, the raw experimental measurements are used without any
preprocessing procedure to attenuate these perturbations. Finally, the wave field contains
waves that are above the breaking limit, as reported by Clauss & Klein (2011), leading
breaking events to be included in the reference data. None of the reconstruction methods
take into account wave breaking and, as explained in the next section, only the third-order
HOS solution (HOSM3) models wave breaking during the surface propagation phase (i.e.
$t > t_0 + T_a$). Hence, the following results also inspect the capability of the reconstruction
approaches to handle wave breaking in the input data, although they do not model it.

5.2. *Adaptation of the algorithms*

Because (i) the observed waves are longer and steeper than in the numerical set-up and (ii)
highly nonlinear potential flow solvers are subject to limitations in terms of wave steepness
and space/modal discretization (Ducrozet et al. 2017), the reconstruction algorithms are
adapted as follows. First, the computational domain of HOS-OCM is extended and
discretized according to its new dimension relative to $\lambda_p$. Considering the limited spatial
extent of the surface measurements, a length $L_x = 2485$ m $\approx 6.8\lambda_p$ is chosen. By keeping
the same number of points per peak wavelength, the grid is made of $N_x = 140$ points.
This extension also affects the choice of the wavenumbers $k_n$ of the wave components
of the analytical models (see § 2.4). Note that the resulting spatial step $\Delta x \approx 17.9$ m is
longer than the wavelength of the measurement perturbations of wavelength $2\Delta x_e = 14$ m
described in the previous section. A similar assimilation period $T_a = 3T_p$ (here $\sim 50$ s) is
used, and the number of time steps $N_t$ is still fixed to 48 per peak period. The location of the
synthetic radar remains similar, i.e. approximately 470 m after the end of the computational
domain. Then, due to the presence of breaking waves, the third-order HOS propagation of
the reconstructed surface leads the surface elevation to be multi-valued if no energy is
dissipated, making the simulation crash. In consequence, the propagation of the HOSM3
solution is performed by an in-house HOS engine that includes a wave breaking model as
described by Wu (2004, § 4.4). It has been ensured, by comparing predicted surfaces for
the steepest case presented in § 4.4.2, that the two employed HOS solvers lead to consistent
results. The other wave models, whose numerical frameworks handle the propagation
of very steep surfaces without crashing, are propagated without modelling any energy
dissipation.
<figure><img src="image_14.png" /><figcaption>Figure 12. Time evolution of the spatial surface profiles, from the reconstruction ($t - (t_0 + T_a) \in [-50, 0]$ s) to the prediction ($t - (t_0 + T_a) > 0$ s) using data type 1. The wave gauge measurements (—, green) are compared with the linear and nonlinear solutions of HOS-OCM, namely HOSM1 (---) and HOSM3 (—). The solutions outside the theoretical prediction zone (······, grey) are coloured in grey. The extreme wave appears at $t - (t_0 + T_a) = 40$ s.</figcaption></figure>
5.3. Results

The application of the algorithms to the experimental data is illustrated through one
specific configuration chosen such that the New Year wave appears in the middle
of the prediction zone 40 s after the surface reconstruction. Similar to the numerical
investigation, the three data types described in § 2.1.1 are used as surface observations.
The calculation of the prediction zone is similar to that described in § 4.1 with the
characteristics of the wave field generated experimentally. Only the regions on the domain
boundaries that are left out of the prediction zone due to the periodicity assumption
(see § 2.5) are kept constant, i.e. approximately equal to 156 m. This choice follows
the observation that the impact of the periodicity assumption on the accuracy of the
reconstruction is very limited, and taking regions spanning one peak wavelength would
be unnecessarily conservative.

The time evolutions of the reference and reconstructed space series of HOSM1 and HOSM3 using data type 1 are presented in figure 12, together with the propagated surfaces up to 100 s of prediction. Both solutions match well the reference surface in the reconstructed domain, that is for $t - (t_0 + T_a) \in [-50, 0]$ s. The distinction between the linear and nonlinear solutions becomes stronger as time increases during the prediction, the linear solution being subject to an increasing phase delay and shape mismatch. In contrast, HOSM3 is able to predict the main features of the surface even after 100 s of prediction. It can also be noted from figure 12 that the definition of the theoretical prediction zone is conservative as the predicted surfaces also correlate, to a certain extent, with the reference outside the boundaries denoted by the dotted lines. A refined procedure for the determination of the actual prediction zone would be of practical interest to avoid underestimating the prediction horizon.

The accuracy is quantified through the calculation of SSP<sub>η</sub> (see § 3.2 for definition) over the reconstructed or predicted domain. In the latter case, the integral terms of (3.1) are bounded by $k_\ell$ and $k_h$ (calculated from $\omega_\ell$ and $\omega_h$, see § 4.1), the limiting wavenumbers of the waves included in the prediction zone. The results are presented in figure 13 for all wave models and data types. In agreement with the numerical results, HOSM3 and ICWM give the most accurate reconstructions/predictions for all data types, with a significant advantage for HOSM3. The linear solutions HOSM1 and LWT follow approximately the same trend and become rapidly incorrect as time increases. The LWT-CDR is consistent with the numerical results in the case of data types 1 and 2, i.e. its reconstruction accuracy is close to that of LWT, but gets closer to that of ICWM during the prediction phase. For data type 3, however, LWT-CDR fails to correctly reconstruct the surface with $SSP_\eta > 0.3$ for all the investigated times. The observed much harder convergence (i.e. need for much more iterations) of the weakly nonlinear AMBM solutions in the case of the experimental set-up compared with the numerical one could explain the depicted inconsistency in the results of LWT-CDR. The determination of the parameter responsible for this lower convergence rate, e.g. the presence of wave breaking events or measurement perturbations, would necessitate further investigations that are outside of the scope of this paper. In contrast, the convergence rate of HOS-OCM is similar in both numerical and experimental configurations.

Finally, figure 14 gives a closer view of the New Year wave prediction for all wave
models and data types. While all the other wave models miss the description of the surface
elevation, HOSM3 is able to predict the correct phase and amplitude of the New Year wave
for all data types. In addition to confirming that all the physics modelled at the third order
is crucial for the correct description of steep wave fields (e.g. Lünser *et al.* 2022), this result
<figure><img src="image_15.png" /><figcaption>Figure 13. Evolution of the spatial SSP of the surface elevation in time during the reconstruction ($t - (t_0 + T_a) \le 0$ s) and prediction ($t - (t_0 + T_a) > 0$ s) periods for HOSM1 (--●--), HOSM3 (—●—), LWT (--Δ-●-, red), LWT-CDR (···Δ···, red) and ICWM (—Δ-—, red). Results using data type (a) 1, (b) 2 and (c) 3 are presented (see table 1 for definitions).</figcaption></figure>

suggests that a prediction system based on such a third-order surface dynamics would be
enough for severe wave conditions.

**6. Conclusions**

The reconstruction and prediction accuracy of surface waves was investigated using
observations representing different modulation mechanisms pertaining to remote
measurements. Using a selection of methods that model different physical properties of
the ocean surface, the impact of the wave shadowing and tilt modulations on the capacity
of the prediction algorithms were highlighted. Reconstructing surfaces from observations
subject to shadowing was found to have an impact limited to the non-visible regions.
However, using radar intensities (assumed proportional to the surface slopes) instead of
surface elevations affected the reconstruction more generally due to the indirect, more
complex extraction of wave information from the measurements. It also brought to light
some numerical properties of the reconstruction methods, HOS-OCM and AMBM, that
influence their respective performances significantly. Still, both methods were able to
predict waves in all configurations with an accuracy that depends on the limitations of the
embedded physical model. It was shown that the accuracy of the surface reconstruction
was mainly related to the correct modelling of the wave shape nonlinearities. In contrast,
modelling the nonlinear correction of the dispersion relation had a marginal effect on the
reconstruction accuracy but substantially helped improve the prediction. The third-order
terms included in HOSM3 that model a frequency-dependent nonlinear wave phase
correction (unlike ICWM) were shown to be crucial for the correct prediction of steep
irregular waves. In addition, the application of the algorithms to data from wave tank
experiments confirmed their suitability for non-perfect measurements, wave breaking and
finite water depth. It demonstrated as well the feasibility of extreme wave predictions with
HOSM3.

The encouraging results of HOS-OCM/HOSM3 call for further investigations on
improving the compatibility of highly nonlinear wave models with reconstruction/
prediction algorithms from real radar measurements and other types of remotely acquired
wave data. Though the presented study is restricted to unidirectional wave fields, the
formalism has been developed to be applicable to directional waves, which will be the
object of future studies. Relying on the explicit resolution of the wave equations also
allows extending of the reconstruction to additional environmental parameters. Ongoing
<figure><img src="image_16.png" /><figcaption>Figure 14. Prediction of the New Year wave, 40 s after the reconstruction. The reference surface (—, green) is compared with the predicted solution obtained using data types 1 (—), 2 (---) and 3 (· · · · · ·) for (a) HOSM1, (b) HOSM3, (c) LWT, (d) LWT-CDR and (e) ICWM.</figcaption></figure>

investigations are considering the inclusion of the current as an unknown of the surface reconstruction problem.

**Acknowledgements.** The authors would like to thank all the reviewers for their constructive comments.

**Funding.** This paper is published as a contribution to the joint research project EproBOSS. The authors wish to express their gratitude to the German Federal Ministry for Economic Affairs and Climate Action (BMWK) and Project Management Jülich (PtJ) for funding and supporting the project (M.K. and N.H., grant no. 03SX510A). Part of this work was funded by the Deutsche Forschungsgemeinschaft (DFG) (N.H., grant no. 277972093).

**Declaration of interests.** The authors report no conflict of interest.
Author ORCIDs.

Nicolas Desmars https://orcid.org/0000-0001-8248-2899;

Moritz Hartmann https://orcid.org/0000-0001-7282-0399;

Norbert Hoffmann https://orcid.org/0000-0003-2074-3170;

Marco Klein https://orcid.org/0000-0003-2867-7534.

Appendix A. Linearization of the radar intensity

Under the assumption that the main property affecting the radar backscatter is the local incidence angle of the radar beam on the surface, it is possible to express the intensity of the backscattered signal by the ocean surface as (Nieto Borge et al. 2004)

$$ \sigma = c_1 \frac{u \cdot n}{|u| |n|} + c_2, \quad (A1) $$

with $c_1$ is a scaling factor and $c_2$ an offset, and where $n$ is the surface-normal vector and $u$ the radar-beam vector, both depicted in § 2.1.1. These vectors are expressed as

$$ n = \begin{bmatrix} -\eta_x \\ -\eta_y \\ 1 \end{bmatrix} = \begin{bmatrix} -\eta_r \cos \beta + \eta_\beta / r \sin \beta \\ -\eta_r \sin \beta - \eta_\beta / r \cos \beta \\ 1 \end{bmatrix}, \quad (A2a) $$

and

$$ u = \begin{bmatrix} x_a - x \\ y_a - y \\ z_a - \eta \end{bmatrix} = \begin{bmatrix} -r \cos \beta \\ -r \sin \beta \\ z_a - \eta \end{bmatrix}, \quad (A2b) $$

from which is deduced

$$ \frac{u \cdot n}{|u| |n|} = \frac{r \eta_r + z_a - \eta}{\sqrt{r^2 + (z_a - \eta)^2} \sqrt{\eta_r^2 + \eta_\beta^2 / r^2 + 1}}. \quad (A3) $$

This leads to a nonlinear relation between $\sigma$ and the surface elevation. By referring as $\epsilon$ to either the surface slopes $\eta_r$ and $\eta_\beta$ or elevation $\eta$, (A3) turns into

$$ \frac{u \cdot n}{|u| |n|} = \frac{r \eta_r}{R} + \frac{z_a}{R} - \frac{r^2}{R^3} \eta + O(\epsilon^2). \quad (A4) $$

As a function of the nominal angle of incidence $\Theta$, it yields

$$ \frac{u \cdot n}{|u| |n|} = \cos \Theta + \eta_r \sin \Theta - \eta \frac{\sin^2 \Theta}{R} + O(\epsilon^2). \quad (A5) $$

The linearized solution of the backscattered radar intensity, assumed valid for the considered sea states, thus reads

$$ \sigma = c_1 \left( \cos \Theta + \eta_r \sin \Theta - \eta \frac{\sin^2 \Theta}{R} \right) + c_2. \quad (A6) $$

Appendix B. Formulation of the HOS-OCM matrices

The formulation of the matrices used in HOS-OCM are described in this section.
## B.1. In the HOS equations

The matrix $\mathbf{A}^{\partial_t}$ is defined according to a finite difference scheme using seven stencils (i.e. of sixth order of accuracy). A central scheme is used, except on the time boundaries of the grid where either forward or backward schemes are used instead. Spatial-derivative matrices $\mathbf{A}^{\partial_x}$, $\mathbf{A}^{\partial_y}$ and $\mathbf{A}^{\partial_z^\ell}$ are derived from the spectral definition of the derivative operators, i.e. $\partial_x \equiv i k_x$, $\partial_y \equiv i k_y$ and $\partial_z^\ell \equiv |k|^\ell$ (with a factor tanh ($|k| h$) for $\ell$ odd), respectively. The notation $\mathbf{A}^{\partial_z}$ in (2.7) refers to $\mathbf{A}^{\partial_z^1}$. The vertical-velocity vector at order $M$ follows:

$$W_M = \sum_{m=1}^{M} W^{(m)}, \qquad (B1)$$

and each order is evaluated as

$$W^{(m)} = \sum_{\ell=0}^{m-1} \frac{\eta^\ell}{\ell!} \mathbf{A}^{\partial_z^{\ell+1}} \cdot \boldsymbol{\phi}^{(m-\ell)}. \qquad (B2)$$

The vector $\boldsymbol{\phi}^{(m)}$ contains the approximation of the velocity potential at $z=0$ at order $m$ and is derived iteratively as

$$\boldsymbol{\phi}^{(1)} = \boldsymbol{\phi}^s, \qquad (B3)$$

$$\boldsymbol{\phi}^{(m)} = - \sum_{\ell=1}^{m-1} \frac{\eta^\ell}{\ell!} \mathbf{A}^{\partial_z^\ell} \cdot \boldsymbol{\phi}^{(m-\ell)}, \quad m = 2, \dots, M. \qquad (B4)$$

To avoid spurious numerical results due to aliasing during the calculation of the products, a full dealiasing procedure is employed (e.g. Bonnefoy et al. 2010). Finally, in order to prevent potential numerical instabilities from appearing in case of high wave steepness, high-frequency modal coefficients of the solution vectors $\boldsymbol{\phi}^s$ and $\eta$ that are used to construct the right-hand side terms $\mathbf{b}^d$ and $\mathbf{b}^k$, and of the right-hand side terms themselves, are damped according to an exponential filter before performing the system inversion.

## B.2. Observation-mapping function

The observation-mapping function relates the quantities at the HOS grid points to the observations. It performs two basic operations, an interpolation (from the grid point locations to the observation locations, in both space and time) and a mapping (from the wave-quantities space to the observation space). More specifically, the latter operation is based on a function that expresses the measured quantity from the surface elevation and/or potential at the grid points. In the presented applications, the observations are related to the surface elevation only, as it is the case for considered wave sensors that give either $\eta$ or $\sigma$ according to (2.2). The interpolation function takes the form of a $N_o$-by-$N_p$ matrix, denoted $\mathbf{A}^i$, that uses Lagrange cubic polynomials. The mapping function depends on the considered observations, and its formulation is detailed below.

### B.2.1. Mapping surface elevations

In the case the observations consist of surface elevations (i.e. data types 1 and 2), the constraint related to the observations is formulated as a simple interpolation following
$η^o = A^i \cdot η$. To match the formalism of (2.7), this equality can be recast into

$$A^M \cdot \begin{bmatrix} \phi^s \\ \eta \end{bmatrix} = \zeta, \qquad (B5)$$

with $A^M = \{O, A^i\}$, $\zeta = \eta^o$, and $O$ an $N_o$-by-$N_p$ matrix that is full of zeros.

### B.2.2. Mapping radar intensities

In the case the observations consist of radar intensities, the constraint related to the observations is formulated based on the linear radar model (2.2). Consequently, the constraint related to the radar observations reads

$$\sigma^o = c_1 \left\{ \cos \Theta^o + \left[ (A^i \cdot A^{\partial_x} \cdot \eta) \cos \beta^o + (A^i \cdot A^{\partial_y} \cdot \eta) \sin \beta^o \right] \sin \Theta^o - (A^i \cdot \eta) \frac{\sin^2 \Theta^o}{R^o} \right\} + c_2 w, \qquad (B6)$$

where vectors having the superscript o contain quantities related to observations, namely the radar intensity $\sigma$, the angle of incidence $\Theta$, the azimuth $\beta$ or the distance $R$ and $w$ is a vector of length $N_o$ that is full of ones. The last equation is equivalent to

$$\frac{\sigma^o}{\sin \Theta^o} = c_1 \left[ \frac{\cos \Theta^o}{\sin \Theta^o} + (D^c \cdot A^i \cdot A^{\partial_x} + D^s \cdot A^i \cdot A^{\partial_y} - D^e \cdot A^i) \cdot \eta \right] + \frac{c_2}{\sin \Theta^o}, \quad (B7)$$

in which $D^c$, $D^s$ and $D^e$ are diagonal matrices built from vectors $\cos \beta^o$, $\sin \beta^o$ and $\sin \Theta^o/R^o$, respectively. To match the formalism of (2.7), equation (B7) is recast into

$$A^M \cdot \begin{bmatrix} \phi^s \\ \eta \end{bmatrix} = \zeta, \qquad (B8)$$

with

$$A^M = \{O, D^c \cdot A^i \cdot A^{\partial_x} + D^s \cdot A^i \cdot A^{\partial_y} - D^e \cdot A^i\}, \qquad (B9)$$

and

$$\zeta = \left( \frac{\sigma^o}{c_1} - \frac{c_2 w}{c_1} - \cos \Theta^o \right) / \sin \Theta^o. \qquad (B10)$$

## Appendix C. Expression of the analytical Jacobian matrices for the least-square problems of AMBM

In this section, the explicit formulation of the Jacobian matrix is given for all configurations of AMBM. Each one is characterized by a wave model (namely LWT, ICWM or LWT-CDR) and an observed quantity (i.e. surface elevation or radar intensity). For the considered wave models, the model parameters are always of the same form, i.e. $p = \{a_1, \dots, a_N, b_1, \dots, b_N\}^T = \{v_a^T, v_b^T\}^T$, where $v_a = a_n$ and $v_b = b_n$ with $n = 1, \dots, N$ and $N$ the number of wave components. The matrix notation is used for conciseness and its similarity to the syntax of matrix-computation programming languages.
### C.1. With elevation observations

Using the surface elevation as observed quantity, the cost function described by (2.11) writes

$$C = \frac{1}{2} |\eta_m(p) - \eta^o|^2, \qquad (C1)$$

with $\eta_m$ the wave model evaluations of $\eta^o$. The corresponding $N_o$-by-$N$ Jacobian matrix is formed by the partial derivatives of every component of $\eta_m$ with respect to every model parameter in $p$. Explicit expressions are detailed below for LWT, ICWM and LWT-CDR.

#### C.1.1. Linear wave theory

The expression of the surface elevation according to LWT is given by (2.13), leading to the Jacobian matrix $J_{\eta,LWT} = \{B^T, C^T\}$ where $B = \cos(\kappa_x \otimes x^o + \kappa_y \otimes y^o - \omega \otimes t^o)$ and $C = \sin(\kappa_x \otimes x^o + \kappa_y \otimes y^o - \omega \otimes t^o)$ with $\otimes$ the outer product, $\kappa_x = k_{x,n}, \kappa_y = k_{y,n}, \omega = \omega_n$, and $x^o, y^o, t^o$ contain the x-, y-, t-coordinates of the observations, respectively.

#### C.1.2. Improved choppy wave model

The expression of the surface elevation according to ICWM is given by (2.15), leading to the Jacobian matrix $J_{\eta,ICWM} = \{E^T, F^T\}$ with

$$E = G - (D^a \cdot H - D^b \cdot G) \odot K + D^{ka} \cdot U, \qquad (C2)$$

$$F = H - (D^a \cdot H - D^b \cdot G) \odot L + D^{kb} \cdot U, \qquad (C3)$$

where

$$K = D^k \cdot [C^{nl} - (D^{ka} \cdot B^{nl} + D^{kb} \cdot C^{nl} + U) \odot (\kappa \omega v_a) \otimes t^o], \qquad (C4)$$

$$L = D^k \cdot [-B^{nl} - (D^{ka} \cdot B^{nl} + D^{kb} \cdot C^{nl} + U) \odot (\kappa \omega v_b) \otimes t^o], \qquad (C5)$$

in which $\odot$ denotes the element-wise matrix multiplication. Here, $D^a, D^b, D^k, D^{ka}$ and $D^{kb}$ are diagonal matrices built from vectors $v_a, v_b, \kappa, \kappa v_a$ and $\kappa v_b$, respectively, where $\kappa = k_n$. $B^{nl} = \cos(\kappa_x \otimes x^o + \kappa_y \otimes y^o - \omega^{nl} \otimes t^o)$ and $C^{nl} = \sin(\kappa_x \otimes x^o + \kappa_y \otimes y^o - \omega^{nl} \otimes t^o)$ are the nonlinear counterparts of $B$ and $C$ that include a correction of the dispersion relation, and $\omega^{nl} = \omega_n^{nl}$. Further, $U$ is an $N$-by-$N_o$ matrix that is full of ones. Finally, $G = \cos Q$ and $H = \sin Q$ with $Q$ a modified phase term defined by

$$ \begin{aligned} Q = \kappa_x \otimes \{x^o - \left[-\left(\frac{\kappa_x v_a}{\kappa}\right)^T \cdot C^{nl} + \left(\frac{\kappa_x v_b}{\kappa}\right)^T \cdot B^{nl}\right] + \kappa_y \otimes \{y^o - \left[-\left(\frac{\kappa_y v_a}{\kappa}\right)^T \cdot C^{nl} + \left(\frac{\kappa_y v_b}{\kappa}\right)^T \cdot B^{nl}\right]\} - \omega^{nl} \otimes t^o. \end{aligned} \qquad (C6) $$
### C.1.3. Linear wave theory with corrected dispersion relation

The expression of the surface elevation according to LWT-CDR is given by (2.18), leading to the Jacobian matrix $J_{η,LC} = \{M^T, N^T\}$ with

$$ M = B^{nl} + (D^{ka} \cdot C^{nl} - D^{kb} \cdot B^{nl}) \odot (\kappa \omega v_a) \otimes t^o, \quad (C7) $$

$$ N = C^{nl} + (D^{ka} \cdot C^{nl} - D^{kb} \cdot B^{nl}) \odot (\kappa \omega v_b) \otimes t^o. \quad (C8) $$

### C.2. With radar intensity observations

Using the radar intensity as observed quantity, the cost function described by (2.11) writes

$$ C = \frac{1}{2} |\sigma_m(\mathbf{p}) - \sigma^o|^2, \quad (C9) $$

with $\sigma_m$ the wave model evaluations of $\sigma^o$. The corresponding $N_o$-by-$N$ Jacobian matrix is formed by the partial derivatives of every component of $\sigma_m$ with respect to every model parameter in $\mathbf{p}$. Explicit expressions are detailed below for LWT, ICWM and LWT-CDR.

#### C.2.1. Linear wave theory

The formulation of the LWT radar intensity is obtained by inserting (2.13) into (2.2), which gives

$$ \sigma_{LWT}(\mathbf{x}, t) = c_1 \sin \Theta \sum_{n=1}^{N} \left\{ -a_n \left[ (k_{xn} \cos \beta + k_{yn} \sin \beta) \sin \psi_n + \frac{\sin \Theta}{R} \cos \psi_n \right] \right. \\ \left. + b_n \left[ (k_{xn} \cos \beta + k_{yn} \sin \beta) \cos \psi_n - \frac{\sin \Theta}{R} \sin \psi_n \right] \right\} \\ + c_1 \cos \Theta + c_2. \quad (C10) $$

The corresponding Jacobian matrix is $J_{\sigma,LWT} = \{R^T, S^T\}$ with

$$ R = -(T \odot C + B \cdot D^e) \cdot D^g, \quad (C11) $$

$$ S = (T \odot B - C \cdot D^e) \cdot D^g, \quad (C12) $$

in which $T = \kappa_x \otimes \cos \beta^o + \kappa_y \otimes \sin \beta^o$, the definitions of $\beta^o, \Theta^o, w$ and $D^e$ follow from § B.2.2, and $D^g$ is a diagonal matrix built from the vector $c_1 \sin \Theta^o$.
### C.2.2. Improved choppy wave model

The formulation of the ICWM radar intensity is obtained by inserting (2.15) into (2.2), which gives

$$ \sigma_{\text{ICWM}}(\mathbf{x}, t) = c_1 \sin \Theta \sum_{n=1}^{N} \left[ \begin{aligned} & -a_n \left( \{[k_{xn}(1-l) - k_{yn}d]\cos \beta + [k_{yn}(1-q) - k_{xn}d]\sin \beta\} \sin \Psi_n \right. \\ & \qquad \left. + \frac{\sin \Theta}{R} \cos \Psi_n \right) \\ & + b_n \left( \{[k_{xn}(1-l) - k_{yn}d]\cos \beta + [k_{yn}(1-q) - k_{xn}d]\sin \beta\} \cos \Psi_n \right. \\ & \qquad \left. - \frac{\sin \Theta}{R} \sin \Psi_n \right) - \frac{\sin \Theta}{2R} (a_n^2 + b_n^2) k_n \end{aligned} \right] + c_1 \cos \Theta + c_2, \quad (\text{C13}) $$

with

$$ \left. \begin{aligned} d &= \sum_{n=1}^{N} \frac{k_{xn} k_{yn}}{k_n} \left( -a_n \cos \psi_n^{nl} - b_n \sin \psi_n^{nl} \right), \\ l &= \sum_{n=1}^{N} \frac{k_{xn}^2}{k_n} \left( -a_n \cos \psi_n^{nl} - b_n \sin \psi_n^{nl} \right), \\ q &= \sum_{n=1}^{N} \frac{k_{yn}^2}{k_n} \left( -a_n \cos \psi_n^{nl} - b_n \sin \psi_n^{nl} \right). \end{aligned} \right\} \quad (\text{C14}) $$

The corresponding Jacobian matrix is $J_{\sigma, \text{ICWM}} = \{V^T, W^T\}$ with

$$ V = \left\{ \mathbf{X} \odot \left[ \mathbf{H} + (\mathbf{D}^a \cdot \mathbf{G} + \mathbf{D}^b \cdot \mathbf{H}) \odot \mathbf{K} \right] - (\mathbf{D}^m \cdot \mathbf{M} \cdot \mathbf{D}^c + \mathbf{D}^n \cdot \mathbf{M} \cdot \mathbf{D}^s) \odot (\mathbf{D}^a \cdot \mathbf{H} - \mathbf{D}^b \cdot \mathbf{G}) - \mathbf{E} \cdot \mathbf{D}^e \} \cdot \mathbf{D}^g, \quad (\text{C15}) $$

$$ W = \left\{ \mathbf{X} \odot \left[ -\mathbf{G} + (\mathbf{D}^a \cdot \mathbf{G} + \mathbf{D}^b \cdot \mathbf{H}) \odot \mathbf{L} \right] - (\mathbf{D}^m \cdot \mathbf{N} \cdot \mathbf{D}^c + \mathbf{D}^n \cdot \mathbf{N} \cdot \mathbf{D}^s) \odot (\mathbf{D}^a \cdot \mathbf{H} - \mathbf{D}^b \cdot \mathbf{G}) - \mathbf{F} \cdot \mathbf{D}^e \} \cdot \mathbf{D}^g, \quad (\text{C16}) $$

in which $\mathbf{X} = -\mathbf{T} + \kappa_x \otimes l \cos \beta^o + \kappa_y \otimes d \cos \beta^o + \kappa_y \otimes q \sin \beta^o + \kappa_x \otimes d \sin \beta^o$. Diagonal matrices $\mathbf{D}^m$ and $\mathbf{D}^n$ are built from vectors $\kappa_x \kappa$ and $\kappa_y \kappa$, respectively, and

$$ d = - \left( \frac{\kappa_x \kappa_y \mathbf{v}_a}{\kappa} \right)^{\text{T}} \cdot \mathbf{B}^{nl} - \left( \frac{\kappa_x \kappa_y \mathbf{v}_b}{\kappa} \right)^{\text{T}} \cdot \mathbf{C}^{nl}, \quad (\text{C17}) $$

$$ l = - \left( \frac{\kappa_x^2 \mathbf{v}_a}{\kappa} \right)^{\text{T}} \cdot \mathbf{B}^{nl} - \left( \frac{\kappa_x^2 \mathbf{v}_b}{\kappa} \right)^{\text{T}} \cdot \mathbf{C}^{nl}, \quad (\text{C18}) $$

$$ q = - \left( \frac{\kappa_y^2 \mathbf{v}_a}{\kappa} \right)^{\text{T}} \cdot \mathbf{B}^{nl} - \left( \frac{\kappa_y^2 \mathbf{v}_b}{\kappa} \right)^{\text{T}} \cdot \mathbf{C}^{nl}. \quad (\text{C19}) $$
### C.2.3. Linear wave theory with corrected dispersion relation

The formulation of the LWT-CDR radar intensity is obtained by inserting (2.18) into (2.2), which gives

$$ \sigma_{LC}(\mathbf{x}, t) = c_1 \sin \Theta \sum_{n=1}^{N} \left\{ -a_n \left[ (k_{xn} \cos \beta + k_{yn} \sin \beta) \sin \psi_n^{nl} + \frac{\sin \Theta}{R} \cos \psi_n^{nl} \right] \right. \\ \left. + b_n \left[ (k_{xn} \cos \beta + k_{yn} \sin \beta) \cos \psi_n^{nl} - \frac{\sin \Theta}{R} \sin \psi_n^{nl} \right] \right\} \\ + c_1 \cos \Theta + c_2. \qquad (C20) $$

The corresponding Jacobian matrix is $J_{\sigma,LC} = \{\mathbf{Y}^T, \mathbf{Z}^T\}$ with

$$ \mathbf{Y} = \{ \mathbf{T} \odot [-\mathbf{C}^{nl} + (\mathbf{D}^{ka} \cdot \mathbf{B}^{nl} + \mathbf{D}^{kb} \cdot \mathbf{C}^{nl}) \odot (\kappa \omega \mathbf{v}_a) \otimes \mathbf{t}^o] - \mathbf{M} \cdot \mathbf{D}^e \} \cdot \mathbf{D}^g, \quad (C21) $$

$$ \mathbf{Z} = \{ \mathbf{T} \odot [\mathbf{B}^{nl} + (\mathbf{D}^{ka} \cdot \mathbf{B}^{nl} + \mathbf{D}^{kb} \cdot \mathbf{C}^{nl}) \odot (\kappa \omega \mathbf{v}_b) \otimes \mathbf{t}^o] - \mathbf{N} \cdot \mathbf{D}^e \} \cdot \mathbf{D}^g. \quad (C22) $$

## Appendix D. Velocity potential of ICWM

In this section, the velocity potential of ICWM in its original Lagrangian framework is derived under deep water assumption. Then, the corresponding surface potential is approximated in the Eulerian framework, leading to the expression that is used for comparison with the reference solution.

### D.1. Velocity potential in the Lagrangian framework

A potential $P$ of a Lagrangian flow is defined such that its derivatives according to a particle instantaneous location $\mathcal{R} = \{x(\bar{x}, \bar{y}, \bar{z}, t), y(\bar{x}, \bar{y}, \bar{z}, t), z(\bar{x}, \bar{y}, \bar{z}, t)\}^T$ (with $\bar{\mathcal{R}} = \{\bar{x}, \bar{y}, \bar{z}\}^T$ the reference particle location defined by the surface at rest) give the particle velocity, i.e. $\nabla_{\bar{\mathcal{R}}}P = \bar{\mathcal{R}}_t$ with $\nabla_{\bar{\mathcal{R}}} = \{\partial_x, \partial_y, \partial_z\}^T$. Because the Lagrangian framework leads to the knowledge of $\bar{\mathcal{R}}_t$ at the particle location $\bar{\mathcal{R}}$ instead of at any arbitrary spatial point, the calculation of $P$ by direct integration of $\bar{\mathcal{R}}_t$ with respect to $\bar{\mathcal{R}}$ is not trivial. Alternatively, demonstrating that a function $P$ is a velocity potential of the Lagrangian flow can be done as follows. By definition, a perfect differential $dP$ can be written in the form

$$ dP = \nabla_{\bar{\mathcal{R}}} P \cdot d\bar{\mathcal{R}} = \nabla_{\mathcal{R}} P \cdot d\mathcal{R}, \qquad (D1) $$

with $\nabla_{\bar{\mathcal{R}}} = \{\partial_{\bar{x}}, \partial_{\bar{y}}, \partial_{\bar{z}}\}^T$. Then, using the following relation

$$ (\mathbf{J} \cdot \mathcal{R}_t) \cdot d\bar{\mathcal{R}} = \mathcal{R}_t \cdot d\mathcal{R}, \qquad (D2) $$

with $\mathbf{J}$ the Jacobian matrix of $\mathcal{R}$ defined by

$$ \mathbf{J} = \begin{bmatrix} x_{\bar{x}} & y_{\bar{x}} & z_{\bar{x}} \\ x_{\bar{y}} & y_{\bar{y}} & z_{\bar{y}} \\ x_{\bar{z}} & y_{\bar{z}} & z_{\bar{z}} \end{bmatrix}, \qquad (D3) $$

it is possible to state that finding a function $P$ such that $dP = (\mathbf{J} \cdot \mathcal{R}_t) \cdot d\bar{\mathcal{R}}$ is a perfect differential (equivalent to showing that $\nabla_{\bar{\mathcal{R}}}P = \mathbf{J} \cdot \mathcal{R}_t$) implies that $\nabla_{\mathcal{R}}P = \mathcal{R}_t$, thus that $P$ is a velocity potential.
The components of the free surface (i.e. at $\bar{z} = 0$) particle displacement according to ICWM is given by Guérin et al. (2019). Noticing that the only differences between the classical second-order Lagrangian solution, such as derived by Nouguier et al. (2015, (5.25)–(5.27)) and Pierson (1961, (8) and (27) with a second-order vertical mean correction), and ICWM are (i) the corrected horizontal particle shift and (ii) the neglected interaction terms between waves of different frequencies, the ICWM components $\mathcal{R}$ are extended here to any water depth $\bar{z}$ to yield

$$ \mathbf{x}(\bar{\mathbf{x}}, \bar{z}, t) = \bar{\mathbf{x}} + \sum_{n=1}^{N} \frac{\mathbf{k}_n}{k_n} (-a_n \sin \tilde{\psi}_n + b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}} + \mathcal{U}_s t, \quad (D4) $$

$$ z(\bar{\mathbf{x}}, \bar{z}, t) = \bar{z} + \sum_{n=1}^{N} (a_n \cos \tilde{\psi}_n + b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} + \frac{1}{2} \sum_{n=1}^{N} (a_n^2 + b_n^2) k_n e^{2k_n \bar{z}}, \quad (D5) $$

where $\tilde{\psi}_n = \mathbf{k}_n \cdot \bar{\mathbf{x}} - \tilde{\omega}_n t$, $\tilde{\omega}_n = \omega_n - \frac{1}{2} \mathbf{k}_n \cdot \mathcal{U}_s$, and $\mathcal{U}_s = \sum_{n=1}^{N} (a_n^2 + b_n^2) \omega_n \mathbf{k}_n e^{2k_n \bar{z}}$, and their time derivatives $\mathcal{R}_t$ follow

$$ \mathbf{x}_t(\bar{\mathbf{x}}, \bar{z}, t) = \sum_{n=1}^{N} \frac{\mathbf{k}_n}{k_n} \tilde{\omega}_n (a_n \cos \tilde{\psi}_n + b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} + \mathcal{U}_s, \quad (D6) $$

$$ z_t(\bar{\mathbf{x}}, \bar{z}, t) = \sum_{n=1}^{N} \tilde{\omega}_n (a_n \sin \tilde{\psi}_n - b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}}. \quad (D7) $$

The corresponding Jacobian matrix, written in the form $\mathbf{J} = \{\mathbf{J}_1, \mathbf{J}_2, \mathbf{J}_3\}$, reads

$$ \mathbf{J}_1 = \begin{bmatrix} x_{\bar{x}} \\ x_{\bar{y}} \\ x_{\bar{z}} \end{bmatrix} = \begin{bmatrix} 1 + \sum_{n=1}^{N} \frac{k_{x_n}^2}{k_n} (-a_n \cos \tilde{\psi}_n - b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} \\ \sum_{n=1}^{N} \frac{k_{x_n} k_{y_n}}{k_n} (-a_n \cos \tilde{\psi}_n - b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} \\ \left( \sum_{n=1}^{N} k_{x_n} (-a_n \sin \tilde{\psi}_n + b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}} + \sum_{n=1}^{N} (a_n^2 + b_n^2) 2k_n \omega_n k_{x_n} e^{2k_n \bar{z}} t + O(\epsilon^3) \right) \end{bmatrix}, \quad (D8) $$
$$
J_2 = \begin{bmatrix} y_{\bar{x}} \\ y_{\bar{y}} \\ y_{\bar{z}} \end{bmatrix} = \begin{bmatrix} \sum_{n=1}^{N} \frac{k_{xn} k_{yn}}{k_n} (-a_n \cos \tilde{\psi}_n - b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} \\ 1 + \sum_{n=1}^{N} \frac{k_{yn}^2}{k_n} (-a_n \cos \tilde{\psi}_n - b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} \\ \left( \sum_{n=1}^{N} k_{yn} (-a_n \sin \tilde{\psi}_n + b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}} \right) \left( 1 + \sum_{n=1}^{N} (a_n^2 + b_n^2) 2k_n \omega_n k_{yn} e^{2k_n \bar{z}} t + O(\epsilon^3) \right) \end{bmatrix}, \quad (D9)
$$

$$
J_3 = \begin{bmatrix} z_{\bar{x}} \\ z_{\bar{y}} \\ z_{\bar{z}} \end{bmatrix} = \begin{bmatrix} \sum_{n=1}^{N} k_{xn} (-a_n \sin \tilde{\psi}_n + b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}} \\ \sum_{n=1}^{N} k_{yn} (-a_n \sin \tilde{\psi}_n + b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}} \\ \left( 1 + \sum_{n=1}^{N} k_n (a_n \cos \tilde{\psi}_n + b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}} \right) \left( 1 + \sum_{n=1}^{N} (a_n^2 + b_n^2) k_n^2 e^{2k_n \bar{z}} + O(\epsilon^3) \right) \end{bmatrix}, \quad (D10)
$$

in which $\epsilon \equiv k\sqrt{a^2 + b^2}$ characterizes the wave steepness. Then, retaining terms up to the
second order in $\epsilon$ and ignoring interaction terms of waves having different frequencies,
such as done for ICWM (Guérin et al. 2019), the components of $J \cdot R_t$ yield

$$
J_1 \cdot \mathcal{R}_t = \sum_{n=1}^{N} \frac{k_{xn}}{k_n} \tilde{\omega}_n (a_n \cos \tilde{\psi}_n + b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}}, \quad (D11)
$$

$$
J_2 \cdot \mathcal{R}_t = \sum_{n=1}^{N} \frac{k_{yn}}{k_n} \tilde{\omega}_n (a_n \cos \tilde{\psi}_n + b_n \sin \tilde{\psi}_n) e^{k_n \bar{z}}, \quad (D12)
$$

$$
J_3 \cdot \mathcal{R}_t = \sum_{n=1}^{N} \tilde{\omega}_n (a_n \sin \tilde{\psi}_n - b_n \cos \tilde{\psi}_n) e^{k_n \bar{z}}. \quad (D13)
$$

A function *P* that satisfies the first-order Lagrangian expansion for an irregular
unidirectional wave field is given by Pierson (1961, (8)) and Nouguier *et al.* (2015,
(3.10)). For the second order, Nouguier *et al.* (2015, (4.46)) give a continuous form of
such a function for a directional irregular wave field, while Pierson (1961, (27)) gives its
discrete formulation for a bichromatic two-dimensional wave field. This latter expression
only involves interaction terms of wave components of different frequencies, which are
neglected in the present formulation of ICWM. Hence, the first-order potential function
that includes the ICWM nonlinear phase correction appears to be a good candidate for the
velocity potential and reads

$$P = \sum_{n=1}^{N} \frac{\tilde{\omega}_n}{k_n} \left( a_n \sin \tilde{\psi}_n - b_n \cos \tilde{\psi}_n \right) e^{k_n \bar{z}}. \quad (D14)$$

Indeed, using the above function, the components of $\nabla_{\tilde{\mathcal{R}}}P$ correspond to those of $\mathbf{J} \cdot \mathcal{R}_t$ (D11) to (D13). Consequently, $(\mathbf{J} \cdot \mathcal{R}_t) \cdot d\tilde{\mathcal{R}}$ is a perfect differential of $P$ to the second order, meaning that $\nabla_{\tilde{\mathcal{R}}}P = \mathcal{R}_t$ and $P$ is a velocity potential of ICWM in the Lagrangian framework.

## D.2. Approximation of the surface potential in Eulerian form

In order to compare the value of the potential with the reference solution at any spatial point, an Eulerian form of (D14) is required. Because the derivation of the complete explicit Eulerian formulation is not trivial and would lead to lose the mathematical simplicity of the Lagrangian formulation, only an approximation the Lagrangian solution is calculated. Similar to the procedure used by Desmars et al. (2020) to evaluate the Eulerian surface elevation of ICWM, the quantities are expressed on the free surface, i.e. at $\bar{z} = 0$, and a change of reference is performed to implicitly take into account the horizontal particle shift $\mathcal{U}_s|_{\bar{z}=0}$. This leads to the following expressions of the surface particle horizontal location and velocity potential

$$\begin{aligned} x(\bar{x} - \mathcal{U}_s|_{\bar{z}=0}t, 0, t) &= \mathcal{X}(\bar{x}, t) = \bar{x} + \sum_{n=1}^{N} \frac{k_n}{k_n} \left( -a_n \sin \psi_n^{nl} + b_n \cos \psi_n^{nl} \right) \\ &= \bar{x} + \mathcal{D}(\bar{x}, t), \end{aligned} \quad (D15)$$

$$P(\bar{x} - \mathcal{U}_s|_{\bar{z}=0}t, 0, t) = \mathcal{P}(\bar{x}, t) = \sum_{n=1}^{N} \frac{\tilde{\omega}_n|_{\bar{z}=0}}{k_n} \left( a_n \sin \psi_n^{nl} - b_n \cos \psi_n^{nl} \right), \quad (D16)$$

where $\psi_n^{nl} = k_n \cdot x - \omega_n^{nl}t$ and $\omega_n^{nl} = \omega_n + \frac{1}{2}k_n \cdot \mathcal{U}_s|_{\bar{z}=0}$. Then, the surface potential of ICWM at any spatial point is evaluated by computing the particle horizontal displacement at its instantaneous rather than its reference location according to

$$\mathcal{P}(\bar{x}) = \mathcal{P}(\mathcal{X} - \mathcal{D}(\bar{x})) \approx \mathcal{P}(\mathcal{X} - \mathcal{D}(\mathcal{X})) = \phi_{\text{ICWM}}^{s}(\mathcal{X}). \quad (D17)$$

Finally, the expression of the surface potential of ICWM in the Eulerian framework is

$$\phi_{\text{ICWM}}^{s}(\mathbf{x}, t) = \sum_{n=1}^{N} \frac{\tilde{\omega}_n|_{\bar{z}=0}}{k_n} (a_n \sin \Psi_n - b_n \cos \Psi_n), \quad (D18)$$

where

$$\Psi_n = k_n \cdot \left[ \mathbf{x} - \sum_{m=1}^{N} \frac{k_m}{k_m} \left( -a_m \sin \psi_m^{nl} + b_m \cos \psi_m^{nl} \right) \right] - \omega_n^{nl}t. \quad (D19)$$

The surface potential of LWT-CDR follows directly by removing the effect of the nonlinear shape correction from (D18), i.e. replacing $\Psi_n$ with $\psi_n^{nl}$.
N. Desmars and others

REFERENCES

AL-ANI, M., CHRISTMAS, J., BELMONT, M.R., DUNCAN, J.M., DUNCAN, J. & FERRIER, B. 2019 Deterministic sea waves prediction using mixed space–time wave radar data. *J. Atmos. Ocean. Technol.* **36** (5), 833–842.

ALPERS, W. & HASSELMANN, K. 1978 The two-frequency microwave technique for measuring ocean-wave spectra from an airplane or satellite. *Boundary-Layer Meteorol.* **13** (1–4), 215–230.

ALPERS, W.R., ROSS, D.B. & RUFENACH, C.L. 1981 On the detectability of ocean surface waves by real and synthetic aperture radar. *J. Geophys. Res.* **86** (C7), 6481–6498.

ARAGH, S. & NWOGU, O. 2008 Variation assimilating of synthetic radar data into a pseudo-spectral wave model. *J. Coast. Res.* **52**, 235–244.

BAKER, A.H., JESSUP, E.R. & MANTEUFFEL, T. 2005 A technique for accelerating the convergence of restarted gmres. *SIAM J. Matrix Anal. Appl.* **26** (4), 962–984.

BELMONT, M.R., HORWOOD, J.M.K., THURLEY, R.W.F. & BAKER, J. 2007 Shallow angle wave profiling lidar. *J. Atmos. Ocean. Technol.* **24** (6), 1150–1156.

BLONDEL, E., BONNEFOY, F. & FERRANT, P. 2010 Deterministic non-linear wave prediction using probe data. *Ocean Engng* **37** (10), 913–926.

BLONDEL-COUPRIE, E., BONNEFOY, F. & FERRANT, P. 2013 Experimental validation of non-linear deterministic prediction schemes for long-crested waves. *Ocean Engng* **58**, 284–292.

BONNEFOY, F., DUCROZET, G., LE TOUZÉ, D. & FERRANT, P. 2010 Time-domain simulation of nonlinear water waves using spectral methods. In *Advances in Numerical Simulation of Nonlinear Water Waves*, pp. 129–164. World Scientific.

CLAUSS, G.F. & KLEIN, M. 2011 The new year wave in a seakeeping basin: generation, propagation, kinematics and dynamics. *Ocean Engng* **38** (14), 1624–1639.

DANKERT, H. & ROSENTHAL, W. 2004 Ocean surface determination from x-band radar-image sequences. *J. Geophys. Res.* **109** (C4).

DESMARS, N., BONNEFOY, F., GRILLI, S.T., DUCROZET, G., PERIGNON, Y., GUÉRIN, C.-A. & FERRANT, P. 2020 Experimental and numerical assessment of deterministic nonlinear ocean waves prediction algorithms using non-uniformly sampled wave gauges. *Ocean Engng* **212**, 107659.

DESMARS, N., HARTMANN, M., BEHRENDT, J., KLEIN, M. & HOFFMANN, N. 2022 Nonlinear reconstruction and prediction of regular waves. In *International Conference on Offshore Mechanics and Arctic Engineering*, Ocean Engineering; Honoring Symposium for Professor Günther F. Clauss on Hydrodynamics and Ocean Engineering, vol. 5B. ASME.

DOMMERMUTH, D. 2000 The initialization of nonlinear waves using an adjustment scheme. *Wave Motion* **32** (4), 307–317.

DOMMERMUTH, D.G. & YUE, D.K.P. 1987 A high-order spectral method for the study of nonlinear gravity waves. *J. Fluid Mech.* **184**, 267–288.

DUCROZET, G., BONNEFOY, F., LE TOUZÉ, D. & FERRANT, P. 2016 Hos-ocean: open-source solver for nonlinear waves in open ocean based on high-order spectral method. *Comput. Phys. Commun.* **203**, 245–254.

DUCROZET, G., BONNEFOY, F. & PERIGNON, Y. 2017 Applicability and limitations of highly non-linear potential flow solvers in the context of water waves. *Ocean Engng* **142**, 233–244.

FISHER, A., THOMSON, J. & SCHWENDEMAN, M. 2021 Rapid deterministic wave prediction using a sparse array of buoys. *Ocean Engng* **228**, 108871.

FUJIMOTO, W. & WASEDA, T. 2020 Ensemble-based variational method for nonlinear inversion of surface gravity waves. *J. Atmos. Ocean. Technol.* **37** (1), 17–31.

GANGESKAR, R. 2000 Wave height derived by texture analysis of x-band radar sea surface images. In *IGARSS 2000. IEEE 2000 International Geoscience and Remote Sensing Symposium. Taking the Pulse of the Planet: The Role of Remote Sensing in Managing the Environment. Proceedings* (Cat. No. 00CH37120), vol. 7, pp. 2952–2959.

GRILLI, S.T., GUÉRIN, C. -A. & GOLDSTEIN, B. 2011 Ocean wave reconstruction algorithms based on spatio-temporal data acquired by a flash lidar camera. In *Proceedings of the Twenty-first International Offshore and Polar Engineering Conference*, pp. 275–282. International Society of Offshore and Polar Engineers.

GUÉRIN, C.-A., DESMARS, N., GRILLI, S.T., DUCROZET, G., PERIGNON, Y. & FERRANT, P. 2019 An improved lagrangian model for the time evolution of nonlinear surface waves. *J. Fluid Mech.* **876**, 527–552.

GUIMARÃES, P.V., ARDHUIN, F., BERGAMASCO, F., LECKLER, F., FILIPOT, J.-F., SHIM, J.-S., DULOV, V. & BENETAZZO, A. 2020 A data set of sea surface stereo images to resolve space-time wave fields. *Sci. Data* **7** (1), 1–12.
HASSELMANN, K., et al. 1973 Measurements of wind-wave growth and swell decay during the joint north sea wave project (JONSWAP). *Dtsch. Hydrogr. Z. A* **8** (12).

HLOPHE, T., WOLGAMOT, H., KURNIAWAN, A., TAYLOR, P.H., ORSZAGHOVA, J. & DRAPER, S. 2021 Fast wave-by-wave prediction of weakly nonlinear unidirectional wave fields. *Appl. Ocean Res.* **112**, 102695.

HUANG, W., LIU, X. & GILL, E. 2017 Ocean wind and wave measurements using x-band marine radar: a comprehensive review. *Remote Sens.* **9** (12), 1261.

HUCHET, M., BABARIT, A., DUCROZET, G., GILLOTEAUX, J.-C. & FERRANT, P. 2021 Nonlinear deterministic sea wave prediction using instantaneous velocity profiles. *Ocean Engng* **220**, 108492.

KABEL, T., GEORGAKIS, C.T. & ROD ZEEBERG, A. 2019 Mapping ocean waves using new lidar equipment. In *Proceedings of the Twenty-ninth (2019) International Ocean and Polar Engineering Conference*, pp. 2558–2562. International Society of Offshore and Polar Engineers.

KLEIN, M., DUDEK, M., CLAUSS, G.F., EHLERS, S., BEHRENDT, J., HOFFMANN, N. & ONORATO, M. 2020 On the deterministic prediction of water waves. *Fluids* **5**(1), 9.

LI, G., WEISS, G., MUELLER, M., TOWNLEY, S. & BELMONT, M.R. 2012 Wave energy converter control by wave prediction and dynamic programming. *Renew. Energy* **48**, 392–403.

LONGUET-HIGGINS, M.S. & PHILLIPS, O.M. 1962 Phase velocity effects in tertiary wave interactions. *J. Fluid Mech.* **12** (3), 333–336.

LUDENO, G. & SERAFINO, F. 2019 Estimation of the significant wave height from marine radar images without external reference. *J. Mar. Sci. Engng* **7** (12), 432.

LUND, B., COLLINS, C.O., GRABER, H.C., TERRILL, E. & HERBERS, T.H.C. 2014 Marine radar ocean wave retrieval's dependency on range and azimuth. *Ocean Dyn.* **64** (7), 999–1018.

LÜNSER, H., HARTMANN, M., DESMARS, N., BEHRENDT, J., HOFFMANN, N. & KLEIN, M. 2022 The influence of characteristic sea state parameters on the accuracy of irregular wave field simulations of different complexity. *Fluids* **7** (7), 243.

LYZENGA, D.R., NWOGU, O.G., BECK, R.F., O'BRIEN, A., JOHNSON, J., DE PAOLO, T. & TERRILL, E. 2015 Real-time estimation of ocean wave fields from marine radar data. In *2015 IEEE International Geoscience and Remote Sensing Symposium*, pp. 3622–3625.

MA, Y., SCLAVOUNOS, P.D., CROSS-WHITER, J. & ARORA, D. 2018 Wave forecast and its application to the optimal control of offshore floating wind turbine for load mitigation. *Renew. Energy* **128**, 163–176.

MEISNER, E., GALVAGNO, M., ANDRADE, D., LIBERZON, D. & STUHLMEIER, R. 2023 Wave-by-wave forecasts in directional seas using nonlinear dispersion corrections. *Phys. Fluids* **35** (6), 062104.

MÉRIGAUD, A. & TONA, P. 2022 A stochastic approach to short-term ocean wave forecasting: preliminary results using data from a remote sensing imaging system. In *International Conference on Offshore Mechanics and Arctic Engineering*, Offshore Technology, vol. 1. ASME.

MORÉ, J.J. 1978 The Levenberg-Marquardt algorithm: implementation and theory. In *Numerical Analysis* (ed. G.A. Watson), pp. 105–116. Springer.

NAAIJEN, P., TRULSEN, K. & BLONDEL-COUPRIE, É. 2014 Limits to the extent of the spatio-temporal domain for deterministic wave prediction. *Intl Shipbuild. Prog.* **61** (3–4), 203–223.

NAAIJEN, P., VAN OOSTEN, K., ROOZEN, K. & VAN'T VEER, R. 2018 Validation of a deterministic wave and ship motion prediction system. In *International Conference on Offshore Mechanics and Arctic Engineering*, Ocean Engineering, vol. 7B, pp. 1–8. ASME.

NAAIJEN, P. & WIJAYA, A.P. 2014 Phase resolved wave prediction from synthetic radar images. In *International Conference on Offshore Mechanics and Arctic Engineering*, Ocean Engineering, vol. 8A. ASME.

NIETO BORGE, J.C., RODRÍGUEZ, G.R., HESSNER, K. & GONZÁLEZ, P.I. 2004 Inversion of marine radar images for surface wave analysis. *J. Atmos. Ocean. Technol.* **21** (8), 1291–1300.

NOCEDAL, J. & WRIGHT, S.J. 1999 *Numerical Optimization*. Springer.

NOUGUIER, F., CHAPRON, B. & GUÉRIN, C.-A. 2015 Second-order Lagrangian description of tri-dimensional gravity wave interactions. *J. Fluid Mech.* **772**, 165–195.

NOUGUIER, F., GRILLI, S.T. & GUÉRIN, C.-A. 2014 Nonlinear ocean wave reconstruction algorithms based on spatiotemporal data acquired by a flash lidar camera. *IEEE Trans. Geosci. Remote Sens.* **52** (3), 1761–1771.

PERLIN, M. & BUSTAMANTE, M.D. 2016 A robust quantitative comparison criterion of two signals based on the sobolev norm of their difference. *J. Engng Maths* **101** (1), 115–124.

PIERSON, W.J. 1961 Models of random seas based on the Lagrangian equations of motion. *Tech. Rep.* New York University.

PLANT, W.J. & FARQUHARSON, G. 2012 Wave shadowing and modulation of microwave backscatter from the ocean. *J. Geophys. Res.* **117** (C8).
QI, Y., WU, G., LIU, Y., KIM, M.-H. & YUE, D.K.P. 2018a Nonlinear phase-resolved reconstruction of irregular water waves. *J. Fluid Mech.* **838**, 544–572.

QI, Y., WU, G., LIU, Y. & YUE, D.K.P. 2018b Predictable zone for phase-resolved reconstruction and forecast of irregular waves. *Wave Motion* **77**, 195–213.

SALCEDO-SANZ, S., NIETO BORGE, J.C., CARRO-CALVO, L., CUADRA, L., HESSNER, K. & ALEXANDRE, E. 2015 Significant wave height estimation using SVR algorithms and shadowing information from simulated and real measured x-band radar images of the sea surface. *Ocean Engng* **101**, 244–253.

SIMANESEW, A., TRULSEN, K., KROGSTAD, H.E. & NIETO BORGE, J.C. 2017 Surface wave predictions in weakly nonlinear directional seas. *Appl. Ocean Res.* **65**, 79–89.

SIMPSON, A., HALLER, M., WALKER, D., LYNETT, P. & HONEGGER, D. 2020 Wave-by-wave forecasting via assimilation of marine radar data. *J. Atmos. Ocean. Technol.* **37** (7), 1269–1288.

STØLE-HENTSCHEL, S., SEEMANN, J., NIETO BORGE, J.C. & TRULSEN, K. 2018 Consistency between sea surface reconstructions from nautical x-band radar doppler and amplitude measurements. *J. Atmos. Ocean. Technol.* **35** (6), 1201–1220.

TRULSEN, K. 2007 Weakly nonlinear sea surface waves – freak waves and deterministic forecasting. In *Geometric Modelling, Numerical Simulation, and Optimization: Applied Mathematics at SINTEF* (ed. G. Hasle, K.-A. Lie & E. Quak), pp. 191–209. Springer.

WANG, G. & PAN, Y. 2021 Phase-resolved ocean wave forecast with ensemble-based data assimilation. *J. Fluid Mech.* **918**, A19.

WANG, G., ZHANG, J., MA, Y., ZHANG, Q., LI, Z. & PAN, Y. 2022 Phase-resolved ocean wave forecast with simultaneous current estimation through data assimilation. *J. Fluid Mech.* **949**, A31.

WEST, B.J., BRUECKNER, K.A., JANDA, R.S., MILDER, D.M. & MILTON, R.L. 1987 A new numerical method for surface hydrodynamics. *J. Geophys. Res.* **92** (C11), 11803–11824.

WU, G. 2004 Direct simulation and deterministic prediction of large-scale nonlinear ocean wave-field. PhD Thesis, Massachusetts Institute of Technology.

WU, J., HAO, X. & SHEN, L. 2022 An improved adjoint-based ocean wave reconstruction and prediction method. *Flow* **2**, E2.

YOON, S., KIM, J. & CHOI, W. 2016 An explicit data assimilation scheme for a nonlinear wave prediction model based on a pseudo-spectral method. *IEEE J. Ocean. Engng* **41** (1), 112–122.

YOUNG, I.R., ROSENTHAL, W. & ZIEMER, F. 1985 A three-dimensional analysis of marine radar images for the determination of ocean wave directionality and surface currents. *J. Geophys. Res.* **90** (C1), 1049–1059.

ZAKHAROV, V.E. 1968 Stability of periodic waves of finite amplitude on the surface of a deep fluid. *J. Appl. Mech. Tech. Phys.* **9** (2), 190–194.