Deterministic wave forecasting with the Zakharov equation

Abstract This study investigates deterministic wave forecasting from the perspective of the Zakharov equation. Forecasts based on linear dispersion, weakly nonlinear amplitude dispersion and the Zakharov equation are compared for reference ocean surfaces generated from the Joint North Sea Wave Observation Project and Pierson–Moskowitz spectra. This approach allows for the role of nonlinearity to be isolated, and demonstrates the success of simple frequency corrections in forecasting wave fields up to moderate steepness. The role of second-order bound waves is investigated by means of analytical formulae, and their impact on forecast accuracy is illustrated for a range of forecast times.


Introduction
Recent years have seen a surge of interest in deterministic wave forecasting for marine decision support systems, a sector that is poised to continue to grow along with further investment in the marine environment. Increasing the efficiency and safety of offshore renewable energy, autonomous shipping and navigation, and supply and offload operations from very large floating structures (see e.g. Zhao et al. 2018) relies on accurate, timely forecasts of wave conditions. Such operations cannot rely on well-developed stochastic wave forecasts, which provide information of a space/time-averaged nature, on typical scales of ∼100 km and ∼6 h. Rather, a deterministic approach is needed, valid over typical scales of ∼1 km and ∼5 min. The information provided about the sea state can then be used to calculate forces and moments on ships or structures (Kusters et al. 2016), inform decisions about the feasibility of marine operations (Belmont et al. 2014), or provide input to adaptive control systems (Fusco & Ringwood 2010;Ma et al. 2018).
Deterministic wave forecasting is a two-step process: (i) the sea-surface elevation (or some proxy thereof) must be measured, and (ii) the spatio-temporal evolution of that measurement must be computed. Both steps present unique theoretical and practical challenges, although herein only those associated with step (ii) will be addressed. Even if precise information about the sea state is available, the fact that waves in deep water are dispersive produces the major limitation to prediction. At lowest order, the phase velocity c p of the waves is given by c p = ω/k for ω 2 = g|k| the deep-water dispersion relation, with g the acceleration of gravity, ω the radian frequency and k the wavenumber. Thus waves of different frequencies propagate at different speeds. Moreover, the energy of a wave group moves at the so-called group velocity c g = dω/dk, and in deep water c g = c p /2. Therefore, waves within a measured area will move out of it at different speeds, and will mix with other waves originating outside that area. The theoretical region that is amenable to prediction is usually called the 'predictable zone'. When the finite steepness of the waves is accounted for, it can be seen that amplitudes also have an influence on the frequencies, altering the dispersion relation, and thereby the predictable zone. All of these factors impact the ability to produce a forecast from measured data.
For usable forecasts, the sea state must be measured remotely, and for most applications the measurement technique of choice is based on marine wave radar (Borge et al. 2004). Many existing applications (e.g. Connell et al. 2015;Hilmer & Thornhill 2016; Al-Ani, Belmont & Christmas 2020) of deterministic forecasting, both research and commercial, make use of X-band radar backscatter (see Young, Rosenthal & Ziemer 1985;Lee et al. 1995;and others), although approaches based on LIDAR have also been advanced (see e.g. Belmont et al. 2007;Grilli, Guérin & Goldstein 2011). X-band radar has the advantage of a large existing hardware base, as well as greater range, reflecting the relatively greater maturity of the technology.
No matter the source, wave measurements inevitably suffer from a number of uncertainties. Wave buoy measurements exhibit deviations from the actual sea state due to the buoy response (McAllister & van den Bremer 2019), while radar is limited by wave shadowing (Plant & Farquharson 2012) (the inability to observe in the radar shadow of a wave crest). While an accurate measurement of the sea-surface height at a single instant, like that furnished by the WASS (wave acquisition stereo system) of Fedele et al. (2013), would be ideal, radar measurements are further impacted by the rotation time of the radar antennae (typically several seconds (Al-Ani et al. 2020)). Moreover, the spacing of the measurement points cannot be determined a priori (Desmars et al. 2020). This issue was recently addressed by Qi et al. (2018b), who provided a very general formulation of the predictable region based on multiple spatio-temporal measurements, both fixed and moving.
For practical forecasts, with a time horizon limited to several minutes, it is imperative that the forecast is produced quickly. For this reason, linear wave propagation, with its low computational demands, has been widely used (Hilmer & Thornhill 2016;Kusters et al. 2016;Al-Ani et al. 2020). This means that such factors as wave-wave interaction, as well as the presence of bound modes, are necessarily missed. Since the importance of these, and the time scale of nonlinear evolution (theoretically −2 times a typical wave period in deep water, where is the wave steepness), depends on the wave steepness, linear theory may be appropriate for short-term predictions of quiescent states for some marine operations (Belmont et al. 2014). Nevertheless, considerable recent research has been devoted to nonlinear deterministic forecasting.
Early work on nonlinear wave forecasting used classical perturbation theory up to second order for unidirectional (Zhang et al. 1996) and short-crested waves (Zhang et al. 1999). This relied on a division of the amplitude spectrum into frequency bands, with assumptions about distinct interactions between wave modes in each band. Subsequent work, starting with Wu (2004), sought to extend the order of nonlinearity and the physics accounted for by using the high-order spectral method (HOS) formulated by West et al. (1987) and Dommermuth & Yue (1987). This method employs Zakharov's (1968) reformulation of the water-wave problem, together with eigenfunction expansion at each order, to account for arbitrary nonlinearity and large numbers of free modes (see e.g. Mei, Stiassnie & Yue 2018).
Since then, numerous studies have appeared on deterministic wave forecasting using HOS, including Köllisch et al. (2018), who supplemented the HOS with a Newton-Krylov method for determining initial conditions, rather than the optimization procedures employed in Wu (2004) or Blondel, Bonnefoy & Ferrant (2010). Qi et al. (2018a) employed HOS to extend the general linear predictable zone of Qi et al. (2018b), as well as comparing with tank data. They emphasize the importance of nonlinearity in reconstructing a sea state, and develop a method to optimize this reconstruction for a given number of wave modes. Klein et al. (2019) subsequently investigated the use of HOS for a variety of different sea states, and used a surface similarity parameter to quantify the accuracy of the forecast. Law et al. (2020) recently used an artificial neural network trained with HOS simulations to produce predictions. Desmars et al. (2020) have developed a nonlinear Lagrangian model (see Guérin et al. 2019) and employed this for forecasting from non-uniformly spaced data from both HOS simulations and experiments.
Other approaches to nonlinear forecasting focus on using particular closed-form equations, such as the narrow-band cubic Schrödinger equation (NLS), and its extensions (Trulsen 2005;Simanesew et al. 2017). Such model equations have the advantages of providing analytical insight and tractability, as well as computational efficiency, albeit with some loss of predictability for broad-banded or short-crested seas. Klein et al. (2020) have compared several nonlinear Schrödinger models to HOS in a variety of conditions, and found that the inclusion of higher-order dispersion in NLS is crucial in achieving accuracy over a variety of conditions.
In the present paper we pursue wave forecasting using the Zakharov equation (Zakharov 1968), a natural extension of the NLS that has not hitherto been used for this purpose. The Zakharov equation includes wave-wave interaction up to third order and yields the NLS and its extensions as narrow-band simplifications. It also serves as the basis for deriving Hasselmann's kinetic equation, which is widely used in (phase-averaged) oceanic wave forecasting models. Owing to a lack of bandwidth restriction and an explicit separation of resonant and non-resonant terms, the Zakharov formulation is simultaneously general and able to furnish detailed insight into the dominant nonlinear interactions for deep-water waves. Moreover, it is possible to distinguish among those nonlinear interactions which give rise to energy exchange, and those which effect only a mutual frequency correction among waves. Isolating this frequency correction yields a computationally simple but remarkably effective forecast. By revisiting the non-resonant contributions, it is likewise possible to isolate the contribution of bound modes to the evolution of a sea state.
In § 2 we describe how a reference ocean surface can be generated computationally to provide the basis for a forecast, and subsequently give a brief overview of some key concepts relating to discrete data points and Fourier transforms, which will be used for forecasting. In § 3 we introduce three possibilities for the time evolution of the forecast. These are linear frequency dispersion, weakly nonlinear amplitude dispersion derived from the Zakharov equation, or weakly nonlinear amplitude dispersion with energy exchange (using the Zakharov equation). The concept of a predictable region, and its dependence on linear and nonlinear forecasting, are discussed in § 3.2. Numerical results are provided in § 4: narrow spectra from the Joint North Sea Wave Observation Project (JONSWAP) are explored in § 4.1, and examples with broad Pierson-Moskowitz (PM) spectra are treated in § 4.2. Then § 4.3 looks at the inclusion of second-order bound modes and their effect on the forecast; and § 4.4 applies the forecasting methodology to reference seas generated by HOS. Finally § 5 provides some concluding remarks and directions for future research.

Synthetic sea surfaces and measurements
2.1. The synthetic sea surface A computational study of the effects of nonlinearity on deterministic wave forecasting requires that realizations of various sea surfaces are available. Samples from these furnish the measurements that drive the forecast, and those forecasts are compared with subsequent spatio-temporal pictures of the sea surfaces. In common with much other work on deterministic forecasting (see e.g. Desmars et al. 2018;Klein et al. 2020;Law et al. 2020), engineering design spectra will be used to generate the synthetic sea surfaces herein. These design spectra provide the basic parameters of the synthetic sea surface. For the time evolution of these sea surfaces, we will largely employ the Zakharov equation, which includes the effects of cubic nonlinearities. To incorporate higher-order nonlinearities, we also use the HOS method developed by Dommermuth & Yue (1987) and West et al. (1987) to propagate our synthetic sea surfaces (see § 4.4).
The following steps towards generating and propagating the sea surface will be outlined in detail below: (i) select a continuous wavenumber spectrum; (ii) discretize the spectrum with a suitably dense number of modes M; (iii) relate the M modal amplitudes to complex amplitudes B i in the Zakharov formulation; (iv) use the B i to write the free surface η = η 1 + η 2 + O( 3 ), including leading-order terms in η 1 and second-order bound modes in η 2 ; and (v) employ the Zakharov equation to calculate the temporal evolution of η above.
The JONSWAP spectrum is widely used in engineering design, and gives a useful representation of energy density for a variety of wind conditions, including for stormy seas. The spectral shape is depending on wavenumber k, with peak wavenumber k p and tunable parameters α, γ and σ . The PM spectrum, a common one-parameter representation of a fully developed sea, is obtained from (2.1) by setting the peak-sharpening parameter γ = 1 and α = 0.0081. It is important to note that any spectrum is essentially a linear concept, drawing on the idea that a free surface can be described as a superposition of sinusoidal wave modes (see e.g. Holthuijsen 2007).
To generate a synthetic sea surface, the wavenumber spectrum Ψ (k) (either JONSWAP or PM) must first be discretized into a large number of wavenumber 'bins' k i , with 2Ψ (k i ) k = a 2 i giving the corresponding mode amplitude, where k is the width of each bin. This dense discretization furnishes the (initial) magnitude of the complex amplitudes |B k | via The phases φ k corresponding to the magnitudes |B k | are arbitrary, and play no role in the energy spectrum Ψ (k), which is predicated on an assumption of a spatially homogeneous and temporally stationary sea. To third order in the wave steepness , the complex amplitudes evolve according to the Zakharov equation where T npqr = T(k n , k p , k q , k r ) are lengthy kernels (Mei et al. 2018), δ qr np is a Kronecker delta function such that and Δ npqr = ω n + ω p − ω q − ω r is a frequency detuning. This equation captures energy exchange and nonlinear frequency correction due to resonant and near-resonant interactions of quartets of waves (see Mei et al. 2018). This is implicit in the derivation, and requires the frequency detuning to be small in an appropriate sense. (For computational implementation of the Zakharov equation, the criterion Δ npqr /ω p < 2 is used, where ω p is the peak frequency of the input spectrum, and an average steepness.) The Zakharov equation thus contains as special cases the third-order Stokes waves, the nonlinear wave-wave interactions of Longuet-Higgins & Phillips (1962), the nonlinear Schrödinger equation (Zakharov 1968) and many other phenomena. The leading-order free-surface elevation η 1 is obtained from the complex amplitudes via where 'c.c.' stands for the complex conjugate of the preceding expression. The next-order surface elevation η 2 is associated with second-order bound harmonics. These arise due to non-resonant interactions, i.e. the linear harmonic wave description gives rise at second order to terms that include sums and differences of wavenumbers and frequencies like in general, such bound modes do not propagate with the usual (linear or nonlinear) dispersion relation. For the same reason, these waves cannot fulfil both wavenumber and frequency resonance conditions. Their contribution to the free surface can be written in terms of complex amplitudes B i and kernels Mei et al. 2018, chap. 4

.A) as
where the shorthand notation Equations (2.5) and (2.7) together yield a free surface with sum and difference harmonics included, corresponding to the O( ) correction in Wehausen & Laitone (1960, equation (27.25)). (For computational implementation, only modes that are in the support of the spectrum are considered, in order to prevent the accumulation of energy in high wavenumbers.) We will omit higher-order bound waves from the description of the free-surface, although these may be calculated from equation (14.3.5) of Mei et al. (2018).
Thus, to obtain a realization of the reference ocean surface with particular statistical properties (such as significant wave height, or spectral peak period), initial magnitudes |B n | are provided by (2.2) and initial phases φ k are chosen randomly from a uniform distribution over [0, 2π). Inserting these initial data into (2.5) yields a leading-order free surface, and inserting into (2.7) allows second-order bound modes to be added if desired. These reference ocean surfaces at t = 0 will be measured (see § 2.2 below), and snapshots at t > 0 (generated from numerical integration of (2.3)) will be compared to forecasts. The above formulation clearly separates the effects of nonlinearity on wave frequency and energy exchange (included in (2.5)) and the effects of nonlinearity on wave shape (included in (2.7)), allowing for the relative importance of bound modes to be examined.

Discrete measurements and the fast Fourier transform
Provided a free surface η(x, t) composed of long-crested unidirectional waves, we take samples at N points corresponding to a set of physical locations x 0 , x 1 , . . . , x N−1 at a specified time. It is simplest to assume that these data points are equally spaced, and so write x n = nL/N, for N samples covering a length L − L/N. Alternatively, for non-uniformly spaced points, a resampling based on interpolation can be used (see Desmars et al. (2020) for an example of how non-uniformly sampled points may be treated). Then the discrete-time Fourier transform (DTFT) of the sequence {y n } N−1 n=0 is defined as The Y( j) are the discrete-time Fourier coefficients, obtained by multiplying the original sequence by the Nth complex roots of unity, and may be efficiently calculated using the fast Fourier transform (FFT).
The inverse transform is given by which yields an N-periodic extension of the original sequence. Moreover, it is possible to construct a continuous signal by replacing nL/N by a continuous variable x: (2.11) For N even, it will be convenient for subsequent computations to rewrite (2.11) in the form (2.12) The first term Y(0)/N is the mean of the sampled points by (2.9) -an approximation to the mean sea surface over the measurement area. The forecast is thus prepared using N/2 − 1 distinct Fourier modes, which should be significantly smaller than, and not a subset of, the M modes in the generated sea.

Deterministic forecasting
3.1. Linear and nonlinear forecasts Section 2.2 describes the process of (a) taking a continuous sea-surface elevation, (b) measuring it at N equidistant points x n = nL/N, (c) computing the FFT of the measured sequence, and (d) using the inverse transform to construct an L-periodic, continuous approximation to the original sea. That approximation will form the basis for a leading-order forecast, with different possibilities for its forward evolution in time. The possibilities to be explored here are: (i) linear dispersion, assuming modes k n have a frequency ω n = √ g|k n |; (ii) Stokes' corrected dispersion, where the linear frequency is corrected by the presence of other waves to Ω n = ω n + p e np T npnp |B p | 2 ; and (iii) third-order nonlinear evolution using the Zakharov equation.
The simplest spatio-temporal forecast is obtained by assuming that propagating modes are governed by the linear dispersion relation ω m = √ g|k m |, which establishes a one-to-one correspondence between wavenumbers k ∈ R + and (positive) frequency. For unidirectional propagation, this allows for a linear forecast ζ L to be constructed directly from the Fourier coefficients via (2.12) to yield (3.1) where the subscript 'L' is used to denote linear dispersion.
Linear theory is predicated on assumptions of small wave steepness . Thus, for steeper waves, a number of nonlinear effects manifest: (i) amplitude-dependent, mutual frequency correction of different wave modes; and (ii) slow energy exchange among resonant and near-resonant wave modes.
The simplest way to exploit the Zakharov formulation is to employ it to calculate an analogue of Stokes' frequency correction, whereby the finite amplitude affects the dispersion relation of waves. This was first demonstrated for a single wave mode by Stokes (1847), for two wave modes by Longuet-Higgins & Phillips (1962), extended to include the effects of surface tension by Hogan, Gruman & Stiassnie (1988), and subsequently to many modes (Stiassnie 1991) and finite depth (Madsen & Fuhrman 2012).
The Zakharov equation (2.3) makes the calculation of the mutual frequency correction algebraically simple: assuming that the magnitude of the complex amplitudes is constant, so |B n (t)| = |B n (0)|, the Zakharov equation can be integrated in time to leading order, as described by Stuhlmeier & Stiassnie (2019). This leads to the compact expression for the nonlinear frequency correction where ω n is the linear frequency, Ω n the nonlinear corrected frequency and For long-crested waves, as those treated herein, the kernel T npnp takes the following simple form: ( 3.4) Comparing (2.12) and (2.5), the rescaling between Fourier amplitudes Y(m) and initial values for the complex amplitudes B m is for N the number of sampling points, and m = 1, . . . , N/2 − 1, so that the measured sea-surface elevation provides all the ingredients for calculating the corrected frequencies.
In turn, this allows for a second forecast, improving upon (3.1): (3.6) Here the subscript 'S' is used to denote Stokes' corrected dispersion. The measured mode amplitudes Y(m) remain constant, but the frequency Ω m , as well as phase and group velocities, depend on those amplitudes through (3.2). Finally, it is possible to include third-order energy exchange between quartets of waves in the forecast, allowing the mode amplitudes to be functions of time. The forecast is then cosmetically similar to (3.1): The subscript 'Z' is used to denote evolution using the Zakharov equation, and the argument of the (complex) terms Y(m, t) contributes the nonlinear frequency correction where Y(i, 0) are the Fourier coefficients obtained from the measurement. This is equivalent to solving a Zakharov equation with N/2 − 1 modes (with N/2 − 1 < M for M the number of modes in the sea). Here the initial amplitudes and phases are obtained from the measurements directly, whereas the initial data for the sea are initialized from a given spectrum. The N surface-elevation measurements give only N/2 − 1 coupled Zakharov equations because the complex amplitudes B i depend on both the free surface and the potential at the free surface (see (2.12a,b) in Stiassnie & Shemer 1984).
3.2. The predictable region Owing to their dispersive nature, waves measured over a given area of ocean surface will travel through -and ultimately out of -that area at different speeds. Thus the group velocities of the longest and shortest waves in the measurement, denoted c g,L and c g,S , respectively, determine the predictable region. Any deterministic forecast is thus time-limited, with the extent of the predictable region determined largely by the region of sea measured, and the group velocities of the waves therein. Because nonlinear effects change the dispersion relation, as encapsulated in (3.2), they also affect the group velocity, as depicted in figure 1. Since the effect of nonlinearity in (3.2) is to increase phase and group velocities versus linear theory (an effect also observed experimentally by Taklo et al. (2015)), in a manner dependent on the amplitudes of other waves in the sea, the extent to which the predictable region varies depends on .
In figure 1 the predictable regions are shown using linear dispersion ( ABC) and nonlinear dispersion ( ABD), with the lines AC, AD and BC, BD calculated for k L and k S , respectively. For low (mean) steepness, the predictable regions are nearly identical, but, as steepness increases, the nonlinear theory allows for a theoretically larger, and slightly shifted, predictable region. (The predictable region increases in size because the nonlinear frequency correction disproportionately affects shorter waves; see the discussion in Stuhlmeier & Stiassnie (2019).) In what follows, we shall compare the forecast and the synthetic 'sea' within the nonlinear predictable region only.

Results
The above sections describe a procedure generating a synthetic sea from a spectrum, sampling that sea at a number of spatial locations, and using those measurements to generate a spatio-temporal forecast. Deriving the measurements from reference ocean surfaces generated from JONSWAP and PM wavenumber spectra as described in § 2.1 allows for the effects of nonlinearity on the forecast to be examined in detail. Moreover, many statistically identical realizations of the ocean surface can be generated to test each forecasting approach. A schematic of our forecasting approach using a JONSWAP spectral shape is given in figure 2. Figure 2(a) shows the JONSWAP energy spectrum Ψ (k). This is used to generate a reference surface or 'sea' η 1 (x, 0), and in figure 2(b) the samples of the reference ocean surface are represented by red circles. As mentioned in § 3.1, it is convenient to choose these measurement points with uniform spacing in the x-direction. In real-world applications, the sampling length and measurement spacing might be constrained by the range and resolution of a shipboard X-band radar, but for our synthetic seas no such restrictions exist. Measuring at N points over a length L means that wave modes k ∈ {2π/L, 4π/L, . . . , 2π(N/2 − 1)/L} can be resolved. Therefore, L must be sufficiently large to resolve the longest waves of interest, and N must be sufficiently large to avoid aliasing errors. The primary effect of selecting L is to determine the size of the x-t domain for which prediction is possible.
Depending on the situation of interest, not all of the sampled modes may be relevantparticularly long or short waves may not have a significant effect on shipping, for example. To this end, a reduced set of modes k ∈ [k L , k S ] is used, which are depicted by the dotted vertical lines in figure 2(c). These represent the longest and shortest waves accounted for in the forecast.
Following § 3.2, the nonlinear predictable region ABC, depicted in figure 2(d), is calculated. For waves travelling in the positive x-direction, measured waves k ∈ [k L , k S ] will propagate through the region bounded by lines AC with slope 1/c g,L and BC with slope 1/c g,S . As time progresses, waves shorter than k S (originating at t = 0 at points x > B) will infiltrate the region ABC from the right of BC, and waves longer than k L (originating at points x < A) will do likewise from the left of AC. In the scenario depicted, at t = 160 s the information in the measurement region AB has completely dispersed, and no prediction can be made. In the examples below, we present the three forecasting techniques of § 3.1 applied to a variety of cases.

Examples based on JONSWAP spectra
In discretizing the spectrum we employ M = 198 wavenumber bins from k = 0.0012 m −1 to k = 0.24 m −1 , with uniform spacing k = 0.0012 m −1 . This gives the reference sea η 1 (without bound waves, calculated using (2.5) as described in § 2.1, and simulated in time with the Zakharov equation ( not a subset of the latter.) As cutoff wavenumbers for the measured waves we take k L = 0.028 m −1 and k S = 0.25 m −1 , so that the longest and shortest waves in the forecast are 224 and 25 m long, respectively (the peak of our JONSWAP spectrum corresponds to 126 m waves). The parameter α can be interpreted as an energy scale, and is proportional to the significant wave height and average slope as seen in table 1 (σ = 0.08, γ = 5 and peak wavenumber k p = 0.05 m −1 are held constant).
In table 1 we use the integral of the spectrum Just as the predictable region depends on the nonlinearity of the sea state, so does the accuracy of the forecast itself. In figure 3 a number of forecasts are depicted compared to a reference sea initialized using a JONSWAP spectrum (as described above). Three mean steepness values = 0.1, 0.15 and 0.2 are shown (see table 1), and the sea at t = 60 s is compared to forecasts using linear dispersion (3.1) and the Zakharov equation (3.7) (figure 3a,c,e), or linear dispersion (3.1) and nonlinear dispersion (3.6) ( figure 3b,d, f ) vertical lines, which are determined from the group velocities of the longest and shortest waves accounted for. For = 0.1, all three forecasts agree reasonably well -there is little difference between linear and nonlinear dispersion, and scant energy exchange within the framework of the Zakharov equation. For = 0.15 the lower phase and group speeds of the linear theory lead to departures from the reference sea. Finally, for = 0.2, the evolving amplitudes become important, with the best agreement seen for the forecast based on the Zakharov equation. Nevertheless, the simpler forecast ζ S continues to display satisfactory agreement.
It is possible to obtain a measure of the quality of the forecast in space and time by considering the difference between reference sea and forecast as in figure 4, which compares the sea and forecasts based on (3.6) and (3.7) throughout the predictable region. While both forecasts initially match the sea (at time t = 0, on the x-axis of both plots), the forecast based on integrating the Zakharov equation remains very accurate until at least t = 60 s, while the forecast based on nonlinear dispersion alone shows significant deviations, comparable to the significant wave height of the sea itself.
The comparisons shown in figures 3 and 4 are each based on a single realization of the reference ocean surface from a JONSWAP spectrum, so that the quality of the forecast will vary depending on those waves measured over the length L at t = 0 s. An aggregate measure of quality is provided by computing the linear correlation between forecast and sea, a simple measure of whether the forecast rises and falls with the reference sea. For two samples x = {x 1 , . . . , x n } and y = {y 1 , . . . , y n }, Pearson's correlation coefficient is written as (4.2) wherex denotes the mean. The correlation coefficient ranges between −1 and +1, with 0 implying no linear correlation. A value of +1 means perfect positive correlation, e.g. crests in one signal coincide with crests in the second, while −1 means perfect negative correlation, where crests in one signal coincide with troughs in the second. Other aggregate measures of fit exist, including the 'surface similarity parameter' developed by Perlin & Bustamante (2016) and used to investigate deterministic wave predictions in Klein et al. (2020).
For the three values of in table 1, the correlation ρ between sea and forecast is computed at three forecasting times t = 30, 60 and 90 s throughout the predictable interval bounded by the nonlinear group velocities (see figure 1). Repeating this procedure with 50 realizations of the sea surface, each initialized with random, uniformly distributed phases, and averaging the correlation, yields the results shown in table 2, whereρ denotes the average correlation. Note that figure 3 is one representative realization for the middle column (t = 60 s) of table 2. To identify high-quality forecasts, those entries with mean x (m)  sea provided by a PM spectrum, which is discretized using 200 wavenumber bins between 0.001 and 0.2. To compare with the results for a narrower spectrum in § 4.1, we select a PM spectrum with significant wave height 5.7 m (see table 1), which yields k p = 0.029 and = 0.056 from (4.1). Since the PM spectrum represents a scenario of wind blowing over unlimited fetch, only one free parameter is available to describe such a fully developed sea. Given the smaller peak wavenumber, the forecasting region must be increased in size in order to sample waves of interest. Consequently L is extended to 2000 m, and N to 300 (hence 149 Fourier modes) to avoid aliasing errors (from L = 1000 m and N = 200 for the narrow spectra in § 4.1). In addition, the smallest and largest wavenumbers of interest for the forecast are chosen as k L = 0.01 m −1 and k S = 0.2 m −1 , so that the longest and shortest waves resolved are 628 m and 31 m in length, respectively. The increase in the measurement region leads to a concurrent increase in the predictable region.
A comparison of the three forecasts developed in § 3.1 for such a PM spectrum is given in figure 5. Owing to the low steepness of the sea (η 1 calculated from (2.5), without bound waves), the linear forecast ζ L proves excellent at t = 60 s -the mean correlation over the predictable region from 50 realizations is 0.92 for the linear forecast, compared to 0.99 for ζ S and 1.00 for ζ Z . For longer times, corrections to the dispersion relation (3.6) become increasingly important -the mean correlation from forecasting 50 realizations of the sea  to time t = 120 s decreases to 0.75 for ζ L , but remains 0.98 for ζ S and 0.99 for ζ Z . For fully developed seas of such low steepness ( = 0.056), the frequency-corrected forecast is nearly identical to the forecast based on the Zakharov equation (3.7).

Bound modes and deterministic forecasting
Thus far only free waves have been considered in simulating the reference sea surfaces and generating forecasts. The main effect of the bound waves is to sharpen wave crests and flatten the wave troughs. An example is shown in figure 6, where a JONSWAP spectrum with = 0.2 (see table 1) has been used to generate free wave modes (denoted η 1 ) and second-order bound modes η 2 as described in § 2.1. Together, η 1 + η 2 give the second-order free-surface elevation.
The presence of bound modes presents a difficulty for forecasting, since part of the measured energy does not propagate with the expected dispersion relation. Figure 6(b) shows Fourier amplitude spectra of all modes η 1 + η 2 , as well as of the free modes η 1 only. Only the latter propagate in accordance with linear or weakly nonlinear dispersion (3.2), depending on the forecasting model used. Nevertheless, the measurement procedure outlined in § 3.1 perforce identifies all of the mode amplitudes (as shown in blue in figure 6b) as free modes. Given an ensemble of well-resolved spatio-temporal measurements, it is possible to construct a wavenumber-frequency spectrum Ψ (k, ω) without assuming the linear dispersion relation. Energy that does not fall on the curved plane ω = √ g|k| in spectral (k, ω) space can then be attributed to bound waves or, for field measurements, ambient currents. For surface elevation measurements at a single time, by contrast, it is theoretically impossible to distinguish free and bound modes. The smallness of the second-order bound mode contributions means that they do not give rise to appreciable errors in forecasting if the wave steepness is low. In a practical setting, a measurement contains free as well as bound modes, like the blue curve η 1 + η 2 in figure 6(a).
The sea state in figure 7 is computed from a JONSWAP spectrum with = 0.2 (see table 1), which is evolved as described in § 2.1, with bound waves added from the analytical expression (2.7). Two forecasts based on measurements at N = 200 points over L = 1000 m are contrasted. Forecast ζ Z (red solid curve) simply applies the procedures of § 3.1 to the measurements, ignoring the fact that bound modes are present in the sea. It thereby misidentifies the bound mode energy as belonging to free modes (see figure 6b). All modes in this forecast ζ Z are evolved using the Zakharov equation forecast (3.7). By contrast, the second forecastζ Z (yellow dashed curve) presents an ideal situation that cannot be realized in practice: it uses a fabricated measurement of the free waves only at t = 0 s, which is evolved to t = 60 s using (3.7) and (3.8), with the analytical bound modes then added from (2.7). It differs from the sea state through its use of fewer modes obtained from the measurements.
In figure 7 it is difficult to distinguish the error incurred by an inability to accurately identify bound mode contributions from that inherent in forecasting a highly nonlinear sea from a limited number of measurements. As noted previously, the forecast quality depends on the measurement area L and the waves therein for a particular realization of the sea surface. Table 3 considers averages over 50 realizations at a variety of forecast times (t = 30, 60 and 90 s) and average wave slopes (JONSWAP spectra with = 0.1, 0.15 and 0.2). It compares the average correlation in the predictable region of x-space between sea and either the forecast ζ Z (where bound modes are misidentified as free), or the forecastζ Z (an ideal scenario where bound and free modes can be distinguished in the measurement at t = 0 s).
The mean correlation for the ideal forecastρ Sea,ζ Z compares favourably with that seen forρ Sea,ζ Z in table 2, where bound modes were not considered. For waves of low   Table 3. Mean correlation over 50 realizations between a sea with bound waves and a forecast that fails to account for bound modes ζ Z , or an artificial forecast that employs analytical bound modes for the measurement ζ Z . The sea is initialized from a JONSWAP spectrum with = 0.1, 0.15 and 0.2 (see table 1), and comparisons are presented at times t = 30, 60 and 90 s.
to moderate steepness with = 0.1 and 0.15, the differences between ζ Z andζ Z are rather small, and the forecasting error introduced by the bound waves is minimal. For steeper waves, however, the differences are somewhat more marked, particularly for long forecasting times.

Comparison with HOS simulations
In order to validate our approach, and study the effects of higher-order nonlinearities, the three forecasting approaches can also be applied to reference seas generated by means of HOS simulations. In this section, we employ the open-source HOS-Ocean solver (Ducrozet et al. 2016) to generate and propagate wave fields in a two-dimensional domain. HOS order M = 5 is chosen to include higher-order nonlinearities beyond those captured by the Zakharov equation, and a JONSWAP spectrum with γ = 5, k p = 0.05 m −1 and H s as in table 1 is used to generate the initial wave field in deep water.
Each initial wave field (at t = 0 s) is sampled at S = 300 uniformly distributed points over a length of L = 2000 m. The initial wave field is subsequently evolved in time with HOS (the details of which may be found in Ducrozet et al. (2007Ducrozet et al. ( , 2016 and Bonnefoy et al. (2010)), while the sampled points are used to generate forecasts using linear dispersion, nonlinear dispersion or the Zakharov equation. One realization of the sea is compared to the three forecasts at t = 60 s in figure 8  steepness ( = 0.1) all forecasts provide good agreement, and up to moderate steepness ( = 0.15) both nonlinear forecasts ζ S and ζ Z perform very well. For yet steeper waves, the mismatch between sea and ζ S and ζ Z is similar in scale to that seen in figure 3. Since each realization of the sea employs uniformly distributed random phases, an aggregate measure of the quality of the three forecasting approaches for HOS-generated sea surfaces (denoted 'Sea') is provided by the mean correlations shown in forecasts marked bold (with mean correlation greater than or equal to 0.8, to indicate a good match) in either table. Comparison with HOS simulations of order M = 3 show no appreciable difference from those with M = 5 presented in table 4, underscoring the relatively small contributions of these higher-order terms on the scales involved.

Conclusions
We have introduced and developed a deterministic wave forecasting methodology based on the Zakharov equation, and used it to investigate the effect of nonlinearity and bound waves on wave forecasts. The Zakharov equation readily yields analytical insight, and the present approach is complementary to recent work employing the higher-order spectral method for reconstruction of nonlinear wave fields (see recent work by Köllisch et al. (2018), Qi et al. (2018a) and Desmars et al. (2020)).
Under an assumption that all interactions take place close to a central wavenumber, the nonlinear Schrödinger equation and its extensions may be derived from the Zakharov equation (2.3) (see Stiassnie 1984). The Zakharov equation thus generalizes the nonlinear Schrödinger models, which have been widely used to study rogue waves (Onorato & Suret 2016), as well as in recent studies of wave forecasting (see e.g. Trulsen 2005;Simanesew et al. 2017;Klein et al. 2020). Because of its lack of bandwidth restrictions, our method is applicable to both narrow JONSWAP spectra (as in § 4.1) as well as broad PM spectra (as in § 4.2).
Nonlinear effects can play an important role in deterministic wave forecasts, the extent of which is dependent on both average steepness as well as forecasting time. For small wave slopes below = 0.1, linear theory is sufficient to provide forecasts for moderately long time -up to approximately 10T p , for T p the peak wave period of the spectrum. Nevertheless, the weakly nonlinear corrections associated with finite-amplitude effectsanalogues of Stokes' frequency correction -can be easily and effectively computed from the sea-surface measurements, and yield a marked improvement over linear forecasts.
For moderately nonlinear wave fields, up to = 0.15, such third-order frequency corrections are likely to be sufficient to provide accurate forecasts, as very little energy exchange between modes can occur within the predictable time frame. For steeper wave fields, the full Zakharov equation should be brought to bear to obtain an accurate forecast. While the theoretical time scale of the Zakharov equation is t 2 ∼ −2 T p , our results indicate that departures from linear theory begin to appear at t 1 ∼ −1 T p , which must be taken into account when preparing a forecast. For the time scales involved in deterministic forecasting, the energy exchange captured by the cubic nonlinearities of the Zakharov equation is sufficient, as demonstrated by comparisons with HOS of order five in § 4.4. The Zakharov formulation yields an analytical expression for the second-and higher-order bound waves. Bound waves are present in the wave field as a consequence of the harmonic description of nonlinear wave motion, and thus do not obey the natural dispersion relation. When sampling a real sea surface at a number of points, the associated Fourier amplitudes will contain such bound mode contributions in addition to free waves. By including appropriate bound modes in the generated reference sea, it is possible to assess the effect of this misidentification of bound wave energy. For the cases considered, with steepness up to = 0.15, second-order bound modes form a relatively minor contribution. Future work will be devoted to accurate identification of bound waves from a wave record, with a view towards furnishing more accurate forecasts for longer times and highly nonlinear seas.
While this work has examined the simplest case of unidirectional propagation, subsequent studies will consider directionally propagating waves from the perspective of the Zakharov equation. In the first instance, this will impact the geometry of the predictable region, thereby altering the linear forecast. As shown by Stuhlmeier & Stiassnie (2019), the extent of the nonlinear frequency correction is also direction-dependent, giving rise to a complex interplay between wave frequency, amplitude and direction in determining the forecast. The inclusion of local wind and a suitable representation of wave breaking -via forcing terms in the Zakharov equation (see Annenkov & Shrira 2009) -would further present an opportunity for comparison with experiments.