OMAE 2018-78367 PHASE-RESOLVED RECONSTRUCTION ALGORITHM AND DETERMINISTIC PREDICTION OF NONLINEAR OCEAN WAVES FROM SPATIO-TEMPORAL OPTICAL MEASUREMENTS

We investigate a nonlinear phase-resolved reconstruction algorithm and models for the deterministic prediction of ocean waves based on a large number of spatio-temporal optical measurements of surface elevations. We consider a single sensor (e.g., LIDAR, stereo-video, etc.) mounted on a fixed offshore structure and remotely measuring fields of free surface elevations. Assuming a uniform distribution of measurement points over the sensor aperture angles, the density of free surface observation points geometrically decreases with the distance from the sensor. Additionally, wave shadowing effects occur, which become more important at small viewing angles (i.e., grazing incidence on the surface). These effects result in observations of surface elevation that are sparsely distributed. Here, based on earlier work by [1], we present and discuss the characteristics of an algorithm, aimed at assimilating such sparse data and able to deterministically reconstruct and propagate ocean surface elevations for their prediction in time and space. This algorithm could assist in the automatic steering and control of a variety of surface vehicles. Specifically, we compare prediction results using linear wave theory and the weakly nonlinear Choppy Wave Model [2, 3], extended here to an “improved” second order formulation. The latter model is based on an efficient Lagrangian formulation of the free surface and was shown to be able to model wave properties that are important to the proper representation of nonlinear free surfaces, namely wave shape and celerity. Synthetic datasets from highly nonlinear High Order Spectral simulations are used as reference oceanic surfaces. Predicted results are analyzed over an area that evolves in time, using the theoretical amount of information assimilated during the reconstruction of the wave field. For typical horizons of prediction, we discuss the capabilities of our assimilation process for each wave model


INTRODUCTION
The availability of real-time phase-resolved wave fields is important to many offshore applications, such as the optimal maneuvering and operations of surface vessels or ocean renewable energy harvesting systems.However, many standard ocean wave models [4,5] only predict phase-averaged quantities based on spectral representations (e.g., directional wave energy spectrum), which do not provide a detailed description of deterministic wave properties (e.g., free surface elevation and fluid pressure fields).Nevertheless, under the linear wave superposition assumption, such spectral representations can be used to generate phase-resolved wave fields, by combining frequency wave component amplitudes with a set of (typically random) phases.Based on this data, the initial wave field can be computed through an inverse Fourier transform (FT) and its propagation in space and time performed using a linear or nonlinear propagation model (e.g., High Order Spectral (HOS) model [6]).
Here, we instead propose a method for directly reconstructing a phase-resolved nonlinear wave field, based on spatiotemporal optical measurements of the ocean surface.Future sea state conditions are then predicted by propagating the reconstructed waves in time and space using a weakly nonlinear model.This problem was studied in some earlier work on the basis of free surface elevation time series measured at fixed wave probes [7][8][9].However, in situ measurements, typically made from a moving vessel or vehicle, are more challenging since wave reconstruction must rely on data acquired at constantly updated locations surrounding the path of the vehicle, which leads to practical limitations.A successful solution to this problem was proposed based on X-band radar measurements made from an onboard sensor, combined with a 3D-FT to reconstruct a large patch of free surface elevation surrounding the sensor [10][11][12][13].This initial estimate was then used in a direct numerical simulation of future sea states.
Our work is based on the new nonlinear ocean surface reconstruction algorithm proposed by Nouguier et al. [1,14], on the basis of simulated free surface elevation datasets, such as would be acquired by an optical sensor that can remotely and simultaneously measure surface elevations at many geo-referenced locations at a high frequency (e.g., a generic LIDAR camera, stereovideo acquisition system).Their algorithm was based on applying the weakly nonlinear Choppy Wave Model [2,3], which is Lagrangian and thus more efficient for propagating wave fields in the time domain than Eulerian models based on higher-order Stokes expansions.In this paper, we both validate and extend this approach by using simulation results from a HOS model as reference data.Linear wave propagation models can be used to provide short-term forecasts in calm sea states.However, as time between the prediction horizon and the initial sea state reconstruction increases and/or the sea state becomes more severe, nonlinear wave effects play an increasingly large role in the proper representation of wave fields, their propagation (e.g., due to am-plitude dispersion effects on wave celerity), and on the accuracy of the forecast.Using a nonlinear wave reconstruction and propagation algorithm thus becomes crucial to an accurate sea state prediction.
In this study, the performance of the weakly nonlinear Choppy Wave Model (CWM) and its improved second-order formulation developed hereby, are compared with predictions of the higher-order, but much more computationally demanding, "HOS-ocean" model [6].Earlier work on ocean surface reconstruction and forecasting algorithms [15][16][17] has shown that the size of the dataset assimilated during the reconstruction phase yields a well-defined time evolving spatial area, over which a sea-state prediction is theoretically possible.Accordingly, in this paper, we also evaluate the efficiency of our prediction method with respect to this theoretical prediction zone.
This paper is organized as follows.We first describe the wave models used for the free surface reconstruction and propagation.Second, for each wave model, we present the inverse problem that must be solved to perform the initial free surface reconstruction.Third, we detail the method for generating synthetic spatio-temporal free surface elevation datasets, such as would be acquired by a the depicted optical sensor.Fourth, for one-dimensional measurements used for simplicity, we analyze the effect on the forecast accuracy of reconstruction parameters, such as high cut-off frequency of the reconstructed spectrum, heterogeneity of the observation grid, and number of spatial observation times.For a two-dimensional surface reconstruction, we finally investigate the effect of directional wave spreading.

WAVE PROPAGATION MODELS
Two types of wave models are applied in this work to reconstruct ocean surfaces, based on a measured spatio-temporal dataset.These are based on: (i) linear wave theory (LWT); and (ii) the first-and second-order Lagrangian CWM, which include some intrinsic nonlinear wave properties.Note that the following developments are done under deep water assumption, but the extension to finite depth is straightforward.

Linear Wave Model
We consider a Cartesian coordinate system (x, y, z), with x and y axes located on the mean water surface and the z axis being vertical and positive upward.A linear ocean surface representation is derived by superposing N individual wave components, propagating in the horizontal plane r r r = (x, y), of amplitude A n , wavelength λ n and direction θ n with respect to the x-axis (n = 1, ..., N), yielding, where t is time, ϕ n are random phases uniformly distributed in [0, 2π], and where k n = 2π/λ n are wavenumbers related to angular frequencies ω n through the deep water dispersion relationship, ω 2 n = k n g, where g denotes the acceleration of gravity.To simplify the mathematical developments, we will use the equivalent formulation, where ψ n = k k k n • r r r − ω n t are spatio-temporal phases, and (a n , b n ) = (A n cos ϕ n , A n sin ϕ n ) are wave parameters describing the ocean surface.

Choppy Wave Model
Optical measurements of the ocean surface provide, at each time, the surface elevation measured for a set of spatial points at irregular but defined Eulerian locations in a reference coordinate system.Hence, the wave model used in the reconstruction algorithm must also be able to provide and use comparable information.On the other hand, Lagrangian representations of the free surface of a given order can model nonlinear wave properties that are not included in Eulerian developments of the same order [2,3,18].Here, the CWM following our developments meets both of these requirements, as it is derived as a solution of Lagrangian dynamical equations in an Eulerian system.
Let R R R L and Z L denote the horizontal and vertical displacements of a particle with initial horizontal position r r r 0 and vertical position z 0 = 0 (i.e., on the free surface), respectively.A Lagrangian perturbation expansion in wave steepness of these parameters yields [19][20][21], where D D D i and Z i are the i th -order terms in the Lagrangian expansion.Expressing the particle coordinates from Eqn. (3) in an Eulerian framework requires finding an explicit relationship between R R R L and Z L , which can be achieved through the implicit relationship [14], The CWM surface elevation based on this expansion are then reformulated as, The zeroth-order Lagrangian solution is the particle position at rest, i.e., D D D 0 = (0, 0) and Z 0 = 0.The first-order solution yields [19][20][21], Assuming a single wave component, this first-order solution has been shown to be consistent with a third-order Stokes expansion for a surface frozen at t = 0 (see Eqn. (63) in [2], with a misprint in the sign of the cos (2Kx) factor).Apart from the mean surface level this makes the first-order CWM efficient at representing nonlinear geometrical properties of ocean surface waves (i.e., sharper crests and flatter troughs).However, at this order, water particles located on the free surface move around circles around their initial position (r r r, z = 0) with a uniform angular velocity and, hence, their mean position is at z = 0. Thus, the asymmetry observed in nonlinear wave shape is not properly accounted, as it would require a non-zero mean surface level [2,18,22].
Increasing the order of the Lagrangian expansion, we find that the second-order term Z 2 includes a vertical shift of the mean water level [19,21], which causes particles to oscillate around a higher position.This important feature was added to the reconstruction process by reformulating the first-order CWM as, η (r r r,t) which, throughout the paper, will be referred to as CWM1.
In addition, we also developed and used an improved formulation for a second-order CWM, denoted ICWM, with, where n is the Stokes drift vector.Note that, for a monochromatic wave, this introduces corrections of the phase and group velocities consistent with third-order Stokes theory.In fact, here, we only consider self-interaction terms at second order, thus neglect the effect of interactions of different components.The detailed derivation of this formulation will be presented in future work [23].

METHODS
We present a theoretical, yet realistic, setup for an optical system able to instantaneously remotely measure a large dataset of ocean surface elevations.This dataset is used in combination with various wave models (e.g., LWT, CWM1, ICWM) to perform a nowcast (i.e., initial) ocean surface reconstruction, and based on this to issue short-term forecasts of ocean surface elevations over a specified area.To this effect, we discuss the suitable choice of cutoff frequencies and directions of the reconstructed wave field, and based on this define the area over which the accuracy of the sea state forecast will be evaluated, as compared to a prediction made with a higher-order HOS model.We detail the free surface reconstruction process for linear and nonlinear wave models, and the regularization procedure that is required for solving the resulting inverse problem.

Setup Description
Let us consider an optical system (denoted by OS throughout this paper) simply mounted on a fixed offshore structure at location (x c , y c , z c ) (z c = 30 m is used in the present applications) and facing the ocean surface with viewing angles (α, β ) = (76 • , 0 • ) and aperture angles (α a , β a ) (here α a = 20 • ), and remotely measuring free surface elevations using J rays (here J = 64 × 64 = 4096), uniformly distributed over the aperture angles.This yields the observation zone shown in Fig. 1.During an assimilation time T a , the OS acquires K sets of spatial observations, at a constant sampling frequency f s (here 1 Hz), leading to J × K = L spatio-temporal observations.
We assume a random incident wave field, with its main di-rection along the x-axis, described by a directional JONSWAP (JS) spectrum, S (k, θ ) = F (k) × G (θ ) [24], of significant wave height H s , peak period T p , and peakedness parameter γ.The directional spreading function is standard and defined as, where θ dir = 0 is the mean direction of propagation of the wave field, and σ is a directional parameter.
The combination of wave field parameters selected in the present simulations is:

Cutoff Frequencies and Directions
In practical applications, the reconstructed wave field components should be limited in frequency and direction to the bandwidths that are meaningful to the dynamic response of the structure or ocean vehicle of interest.As the study of a specific structural response is beyond the scope of this paper, we first aim at accurately reconstructing wave fields without specifying any response bandwidth.Hence, the cutoff limits k min,max , θ min,max are calculated with respect to the reference ocean wave spectrum.
Additionally, both the size of the observation area and the spatial sampling resolution are related to the minimum and maximum wave frequencies that can be reconstructed, respectively.For the sake of simplicity, here, we only consider spatial observations (i.e., K = 1) to estimate cutoff frequencies.In this case, the low cutoff wavenumber (or frequency since ω = √ gk) is defined by the largest distance L obs between two observation points as, k min ≥ 2π/L obs .In the following, we set k min to its minimal value by determining L obs at the first observation time t obs .When reconstructing a signal over a regular observation grid (i.e., with constant spatial sampling frequency), the maximum high cutoff frequency must satisfy Shannon's condition, i.e., k max ≤ 2π/ (2 obs ), where obs is the distance between two observation points.Since the observation grid is highly irregular, we will determine k max by evaluating its influence on the quality of the reconstruction.
Cutoff limits in wave directions are set by considering the fraction of total wave energy beyond which the remaining energy can be considered as negligible, i.e., with µ the fraction of total energy considered as negligible.In the following µ = 1%, which corresponds to θ min,max ±35 • for the spreading function defined above.

Prediction Zone Definition
Only a limited amount of information is assimilated during the reconstruction process, due to the physical spatio-temporal limitations of optical free surface observations (e.g., limited size of observation zone) and the finite bandwidths of the reconstructed wave field.Hence, only a similarly limited spatiotemporal region, called prediction zone, is accessible for sea state forecast, as a result of propagating the assimilated information [16,17].For a one-dimensional (1D) wave field, this region is bounded by the time evolution of the amplitude and phase of the slowest and fastest wave components at the end (x max obs ) and beginning (x min obs ) of the observation zone, respectively.It can be shown that reconstructed information attached to a specific wave component (i.e., wave amplitude and phase) is traveling at the wave group velocity C g = ω/ (2k) [15].Thus, in 1D, as shown in Fig. 2a, a point (x,t ≥ t rec ) is included in the prediction zone if, where t rec and t obs are the reconstruction and observation times, respectively.Both the LWT and CWM1 models rely on the linear deep water dispersion relationship, while ICWM includes a second order correction term which affects both phase and group velocity.The two constrains described by Eqn. ( 12) eventually cross each other (at point P max = (x max ,t max ); Fig. 2a), which means that the measured information is obsolete and can no longer be used in a prediction.Note that increasing the assimilation time T a leads to increasing the size of the prediction zone.
In two-dimensions (2D), the wave field is decomposed in individual wave components propagating in direction θ ∈ θ min , θ max .The corresponding prediction zone can be calculated for each direction using Eqn.(12), as if it were a 1D case, by replacing x obs by d obs = r r r obs • k k k = x obs cos θ + y obs sin θ , the distance along the considered direction.The intersection of each zone forms the 2D prediction zone.To save computational time, the prediction zone area can be approximated based on the two extreme directions θ min and θ max , as shown in Fig. 2b.The farthest reachable prediction zone P max corresponds to a segment (Fig. 2b).

Data Assimilation Process
Linear Reconstruction For a linear wave field represented by Eqn.(1), the reconstruction process is based on the minimization of a quadratic cost function representing the mean square difference between observations and predictions of the wave model.Here, assuming there are no measuring and model errors, the cost function is, where η (r r r ,t ) and η are the predicted and measured surface elevations at observation points r r r and at time t ( = 1, ..., L), respectively.This problem is equivalent to a least-square minimization of a cost function, which can be achieved by specifying (for n, m = 1, ..., N; note that index summation is implied for repeated indices), where, is the control vector of 2N unknown wave parameters and, is a vector containing the observation points information, with

and,
is the wave model matrix.The linear system of Eqns.( 14) to (17) can then be solved for the wave model parameters (a n , b n ).
Nonlinear Reconstruction For a nonlinear sea state reconstruction the cost function of Eqn. ( 13) is also minimized, using one of the nonlinear wave models, CWM1 or ICWM, to represent surface elevations (Eqns.(7) or (8)).Here, we detail the system construction for CWM1.Following the same methodology as for the linear model, we use the cost function minimization Eqn. ( 14) with, with, Similarly, we have, Since both A mn and B m now depend on wave parameters (a n , b n ), the system of Eqns.( 14) becomes nonlinear and is solved iteratively.To do so, equations are linearized by computing A  16) and (17).Convergence is typically achieved within 5 to 10 iterations.
A similar approach can be used for ICWM, but the inverse problem becomes more complex due to the higher-order of expansion in the formulation.The complete derivation of the solution for this model and its reconstruction system will be presented in future work.

Regularization
In applications, the inverse reconstruction problem defined above (Eqn.( 14)) can become illconditioned due to practical constraints, such as the heterogeneous distribution of spatial observation points, the limited ocean area observed by the OS, and the the frequency and direction bandwidth cutoffs in the reconstructed wave field.To be able to find consistent results independently of the conditioning of the matrix to invert (i.e., A mn ), a Tikhonov regularization method is applied, which yields, where r is the regularization parameter and ||.|| denotes the Euclidean norm.The optimal regularization parameter is found using the "L-curve" method, which consists of finding the r value corresponding to the point of maximum curvature (i.e., corner) of the (log ||A mn X n − B m || , log ||X m ||) function, hence, providing an optimal compromise between minimizing the residual error and ensuring that the norm of the solution does not become too large.The L-curve corner can be determined analytically through solving a singular value decomposition problem [25,26].

RESULTS AND DISCUSSION
In the following applications, we illustrate and discuss key characteristics of our proposed ocean wave reconstruction algorithms, for optically measured surface elevation data.In the absence of field data, the algorithms are applied to reference synthetic surface elevation datasets, generated first in 1D using a higher-order HOS model.The influence of reconstruction parameters on the accuracy of the solution with respect to the reference data is studied, in particular, the number of wave components, the high cutoff frequency, the heterogeneity of the observation grid, and the number of observation times, while highlighting differences between linear and nonlinear prediction results.We then present a 2D case and show how the horizontal aperture angle can help in the reconstruction.

Synthetic Dataset and Error Definition
The open source HOS-ocean model [6] is used to generate synthetic reference ocean surfaces from which datasets of observation points are extracted.Since its first developments [27,28], this model, which simulates the nonlinear propagation of directional wave fields over large spatio-temporal domains, has undergone numerous validations and been used in various applications, such as for simulating freak waves [29], developing a numerical wave tank [30], or in wave prediction [31,32].Based on an initialization with a JS spectrum and a set of random phases, HOSocean computes the time evolution of the complex amplitude of wave components.These are then used to construct snapshots of the reference ocean surface at selected observation times.In the presented case, the order of the HOS expansion is set to 1 or 5, to generate a linear or a nonlinear reference surface, respectively.In each case, the HOS numerical grid and parameters are selected to ensure converged results.
Based on the computed reference surfaces, synthetic datasets of surface elevations data such as would be acquired by the depicted OS are generated by computing the intersection between J optical rays with the reference ocean surface (see, e.g., Fig. 3 for a 1D surface).
The following relative standard error (RSE) is used as an indicator of the prediction accuracy, where r r r defines the footprint of the spatial prediction zone at time t and N r is the number of reconstructed points, which is selected large enough to achieved a converged error estimate.RSE quantifies the scaled mean squared difference between the reconstructed and reference surfaces over the observation points in the prediction zone.Since the quality of the prediction depends on the local surface geometry at the observation time (e.g., through shadowing effects), an overall RSE is computed by averaging RSE values for N s = 50 surface reconstructions of ocean surface generated at the same time with HOS-ocean for different initial random phases, using the same initial wave spectrum.

One-Dimensional Reconstruction
Using a JS spectrum with parameters as defined above, 1D linear surfaces are first generated and then reconstructed based on data acquired at a single time (i.e., K = 1), with L = 64 observation points uniformly distributed over the observation area.Figure 4 shows the RSE at reconstruction time t rec = t obs = 0 (i.e., for a nowcast), as a function of the number of wave components N used in the reconstructed wave field, for different high cutoff frequencies.Note that Shannon's criterion yields, k max ≤ 13.9k p , but the HOS surface was generated using a cutoff frequency k max 16k p .As expected, results show that the reconstruction error decreases as the high cutoff frequency grows, and eventually converges to a minimum, here approximately 10%, for k max = 10k p .In all cases, the number of wave components (N = 15 to 60) appears to have a very limited influence on the accuracy of the reconstruction.However, this number must be

FIGURE 5. PREDICTION ERROR OF A NONLINEAR OCEAN SURFACE AS A FUNCTION OF TIME, FOR DIFFERENT WAVE RECONSTRUCTION MODELS, OVER A REGULAR GRID FOR
large enough to allow for a proper inversion of the reconstruction system.We find: N ≥ 30 for k max ≥ 8k p and N ≥ 40 for k max ≥ 12k p .In the following we use k max = 10k p and N = 50.
Next, nonlinear reference surfaces are generated using the same parameters and the reconstruction error over the prediction zone is computed as a function of time t = t rec = 0 to 3T p 0.8t max , for the same L = 64 uniformly distributed observation points, using either the LWT, CWM1, or ICWM model  in the reconstruction.For each model, the ocean surface is first reconstructed at t = t rec = 0 and then propagated to a later time where its RSE error with respect to the reference HOS solution is evaluated; for comparison, the case of Fig. 4 propagated in time is also marked on the figure (i.e., reconstruction using LWT of a linear reference surface).Figure 5 shows that at t = t rec = 0 all models give the same results, which is expected since at initial time the different reconstruction algorithms all solve equivalent inverse problems.Then, as time increases, each wave model behaves according to its own properties.While the CWM1 model provides a better prediction than the LWT model, the ICWM further improves the prediction accuracy as time increases, due to the additional nonlinear phase speed corrections.Here, the improvement of the nonlinear over the linear prediction is on the order of a few percent.The linear prediction of a linear reference surface, shown on the same figure, is an indication of the best reconstruction accuracy that could be achieved for this problem and dataset, if all nonlinear effects were properly accounted for in the nonlinear wave models.Hence, the difference between this linear and nonlinear forecasts using ICWM increases with time due to the incomplete representation of time-dependent nonlinear effects in the latter model, to approximately reach 5% at t = 3T p .
Ocean observations made by an OS such as a LIDAR camera yield highly irregular grids, since the density of observation points geometrically decreases with the distance from the OS, which is made worse by wave shadowing effects becoming more important at grazing incidence angles (i.e., for the most distant observation points; see Fig. 3).The effect of using an irregular grid on the accuracy of the ocean surface reconstruction is assessed by recalculating the case of Fig. 5, using the ICWM model, for a dataset of L = 64 irregularly spaced observation points.First, considering a single observation time (i.e., K = 1), Fig. 6 shows that the reconstruction error greatly increases (by more than a factor of 2) as compared to a regular dataset.As suggested in earlier work [1], however, this deterioration of the ocean reconstruction and forecasting can be significantly reduced by increasing the number of observation times.Indeed, this has the effect of making some portions of the free surface, invisible to the OS at a given time due to shadowing effects, visible and hence measurable at a later time.Thus, for K = 2, Fig. 6 shows that the error already drops by 10% as compared to K = 1, and then further drops when increasing K to reach a maximum of ∼ 20% for K = 10.The latter, however, is still 5% larger than the error achieved when using a regular grid, for the same number of observation times (also shown in the figure).Comparing results in Figs. 5 and 6 we also see that, when using a regular grid of observation points, increasing K has no significant effect on the prediction error; for both the LWT and ICWM models, maximum errors indeed remain on the order of 15 to 20% for K = 1 and 10, even if ICWM leads to a consistently better prediction.Figure 6 also shows that the reconstruction error decreases when K increases form 1 to 10, but then gradually increases for larger K values.This might result from increasing numerical errors in the inverse problem solution, as the observed dataset becomes larger.Moreover, since reconstructed waves are only modeled up to second-order in nonlinearity in ICWM, higherorder time-dependent nonlinear effects included in the HOSocean solution, which affect the reconstruction process when using temporal data, are not well represented in the wave model as time increases.Therefore, the model becomes unable to properly represent the ocean surface when the assimilation time T a is too long (e.g., phase shifts resulting from amplitude dispersion effects).Varying the temporal sampling frequency (set here to 1 Hz) would allow finding out whether it is the assimilation time or the total number of observations L that is more likely to improve, and then limit, the quality of the prediction.
In general, increasing the assimilation time T a results both in a decrease of the prediction error and in an increase of the size of the prediction area.Both effects are illustrated in Fig. 7, which shows the reference and reconstructed surface elevations at reconstruction time t rec , for K = 1 and 30.Observation points are shown at the first observation time t obs .With K = 30, the reconstructed surface is much closer to the reference data than for K = 1, and the prediction area is larger.The figure shows that, for K = 30, the reconstructed surface remains quite close to the reference data even beyond the prediction zone (x/λ p ≥ 3.8).This suggests that the estimate for the upper limit of the prediction zone, which is calculated based on the distance traveled by the slowest wave component during the assimilation time, is conservative (Eqn.( 12)).Since longer waves travel faster, the prediction of low frequency components can still be accurate beyond the upper limit.A similar effect affects the high frequency components at distances below the lower limit, when time increases (t > t rec ).This suggests the presence of a "gray area" surrounding the prediction zone and growing with time, over which the amount of measured information could still be enough for an accurate prediction, but with decreasing accuracy as distance increases from the prediction zone boundary.
In Fig. 7, when the density of observation points decreases (K = 1; 1.2 ≤ x/λ p ≤ 2.3), the reconstructed surface tends not to fit well the reference surface elevation.As Shannon's condition is not met over this part of the observation zone, an aliasing phenomenon leads to overestimating high frequency component amplitudes, which causes the observed high frequency surface oscillations.However, since this poorly reconstructed part of the observation zone is located close to the lower limit of the prediction zone, and the inaccuracy only affects high frequency components, the growing "gray area" effect tends to remove that part of the surface from the prediction zone, while the other part remains within it.This explains the decreasing prediction error as a function of time for small K values observed in Fig. 6.With several observation times, this phenomenon disappears due to an artificial refinement of the observation grid: as indicated before, parts of the surface that were initially not illuminated by the OS rays become visible at a different observation time.

Two-Dimensional Reconstruction
We now generate a 2D nonlinear wave field with HOSocean, for the same input JS spectrum as before, and reconstruct it using an observation grid of 64 × 64 spatial points (L = 4096) with k max = 10k p , K = 10, and 50 wave frequencies distributed over 30 directions, yielding N = 1500 wave components.Since the number of spatial observations is significantly larger than before, calculating the exact intersection of L optical rays with time-varying reference surfaces becomes numerically demanding.Thus, we estimate the horizontal location of observation points, r r r , as if the reference surface were plane, neglecting shadowing effects.Results are comparable since the major contribution to the spatial spreading of observation points, i.e., the increasingly grazing incidence of OS rays at large distances, is accounted for.Also, the reconstruction error is now calculated for one surface and not averaged over 50 surfaces as in 1D, and the study presented here is limited to the linear wave reconstruction, even if the reference oceanic surface is nonlinear.
In the following, we estimate the effect of the directional spreading of the observation grid on the ocean reconstruction, which is function of the horizontal aperture angle β a of the OS (Fig. 1). Figure 8 shows the reconstruction error calculated as a function of time over the 2D prediction zone.With β a = 10 • , the error is ∼ 70% at reconstruction time and almost monotonically decreases as time increases up to t = 3T p .Similar to 1D cases, this decrease in error results from the growing "gray area" surrounding the prediction zone boundary.Poorly reconstructed parts of the ocean surface close to the initial prediction zone boundary are gradually removed as time increases.However, with the decreasing size of the prediction zone, the error estimate becomes more subject to slight variations, depending on which waves are crossing the prediction zone at the considered time; this explains the increasing error at t = 2.5T p .Such effects did not occur in the 1D cases since errors were averaged over 50 different surfaces.Figure 8 shows that increasing β a results in a significantly reduced error, which converges to a minimum on the order of 15% for β a = 110 • ; this error is similar to 1D cases.This large error reduction is related to the area of the ocean surface covered by the observation zone, as compared to the area of the prediction zone, at reconstruction time.As shown in Fig. 2b, the two-dimensional prediction zone depends on the range of directions of propagation of wave compo-nents in the reconstructed wave field.Here, as θ min,max ±35 • and the OS is facing the main wave direction, the optimal aperture angle that leads to a maximal prediction zone coverage is 180 − |θ max | + θ min 110 • .Figure 9 shows the horizontal distribution of the local error |η − η HOS | /H s , where η is the reconstructed linear surface and η HOS the nonlinear reference surface, at reconstruction time and after two peak periods of propagation, using an aperture angle β a = 110 • .The prediction error is clearly much smaller within the prediction zone than outside of it.At both selected times, the "gray area" with reduced local error, which surrounds the prediction zone boundary, can clearly be identified.At reconstruction time, the "gray area" is close to the the upper limit of the prediction zone (right side of the boundary in Fig. 9 top) because the reconstructed longer waves propagate faster, and thus reach beyond this limit during assimilation time.After 2T p , the "gray area" has expanded and is now completely surrounding the prediction zone.This confirms that the method used here for estimating the prediction zone boundary is conservative.

CONCLUDING REMARKS
We showed that both 1D and 2D ocean wave fields (or sea states) can be accurately reconstructed on the basis of irregular datasets of realistic, but synthetic (and thus error free), optical measurements.Reconstruction algorithms based on linear and various nonlinear wave models were applied and validated.Despite its simple analytical formulation, the Choppy Wave Model (CWM) is able to simulate relevant nonlinear effects in ocean wave fields.The new improved Choppy models that better predict wave shape (CWM1) and nonlinear phase velocity (ICWM), allow improving the accuracy of the wave reconstruction and prediction.
We quantified effects on the reconstruction of the main characteristic of optical measurements, which is to yield highly sparsely distributed datasets, making for a particularly challenging spatial reconstruction.However, we showed that using spatio-temporal information allows overcoming this issue.
We quantified the influence of several reconstruction parameters, which are shown to be of significant importance to the prediction error.In particular, adequately selecting the cutoff limits of the reconstructed spectrum (in both frequency and direction) and the optical system viewing angles, allows accurately predicting the wave field up to several peak periods of propagation.The influence on the prediction results of a growing "gray area" surrounding the prediction zone boundary was highlighted.
The numerical algorithms used in the present work were not optimized for computational efficiency nor parallelized.On the basis of theoretical and practical estimations, it is expected that optimized and efficiently parallelized algorithms (e.g., on multicore or GPU processors) will easily cope with the increased complexity and computational demand resulting from including non-linear effects in a 2D ocean reconstruction.At present, however, the assessment of such computational performance of the method is beyond the scope of this paper and will be left out for future work.

FIGURE 1 .
FIGURE 1. GEOMETRICAL PARAMETERS LEADING TO THE OBSERVATION ZONE.

FIGURE 2 .
FIGURE 2. (a) SPATIO-TEMPORAL EVOLUTION OF THE PRE-DICTION ZONE FOR A 1D SET OF OBSERVATIONS, AND (b) LO-CATION OF THE PREDICTION ZONE FOR A 2D SET OF OBSER-VATIONS FOR A SINGLE TIME RECONSTRUCTION.

m
based on wave parameters obtained at the previous iteration p, when solving for X (p+1) n at iteration p + 1.The solution is initialized at p = 1, using A (0) mn and B (0) m based on the linear reconstruction problem and thus on Eqns.(

FIGURE 3 .
FIGURE 3. 1D SPATIAL SAMPLING OF A NONLINEAR FREE SURFACE.STRAIGHT LINES ARE OS RAY TRAJECTORIES.

FIGURE 4 .
FIGURE 4. LINEAR RECONSTRUCTION ERROR OF A LINEAR SURFACE AT TIME t rec = t obs = 0 (NOWCAST), AS A FUNCTION OF THE NUMBER OF WAVE COMPONENTS, FOR DIFFERENT HIGH CUTOFF FREQUENCIES, OVER A REGULAR GRID FOR K = 1.

FIGURE 7 .
FIGURE 7. ELEVATION AT RECONSTRUCTION TIME OF THE REFERENCE (HOS) AND RECONSTRUCTED (ICWM) SURFACES, FOR DIFFERENT NUMBERS OF OBSERVATION TIMES.VERTICAL LINES DEFINE THE SPATIAL LIMITS OF THE INITIAL PREDICTION ZONE.

FIGURE 8 .
FIGURE 8. PREDICTION ERROR OF A NONLINEAR OCEAN SURFACE AS A FUNCTION OF TIME, USING ICWM FOR DIF-FERENT HORIZONTAL APERTURE ANGLES.

FIGURE 9 .
FIGURE 9. LOCAL 2D LINEAR PREDICTION ERROR |η − η HOS | /H s , AT RECONSTRUCTION TIME (TOP) AND AFTER 2T p OF PROPAGATION (BOTTOM).LINES MARK THE BOUNDARY OF THE PREDICTION ZONE.
which corresponds to a realistic moderately steep sea state, frequently observable in open oceans, with a corresponding, deep water characteristic steepness, H s /λ p 2%.