An improved adjoint-based ocean wave reconstruction and prediction method

Abstract We propose an improved adjoint-based method for the reconstruction and prediction of the nonlinear wave field from coarse-resolution measurement data. We adopt the data assimilation framework using an adjoint equation to search for the optimal initial wave field to match the wave field simulation result at later times with the given measurement data. Compared with the conventional approach where the optimised initial surface elevation and velocity potential are independent of each other, our method features an additional constraint to dynamically connect these two control variables based on the dispersion relation of waves. The performance of our new method and the conventional method is assessed with the nonlinear wave data generated from phase-resolved nonlinear wave simulations using the high-order spectral method. We consider a variety of wave steepness and noise levels for the nonlinear irregular waves. It is found that the conventional method tends to overestimate the surface elevation in the high-frequency region and underestimate the velocity potential. In comparison, our new method shows significantly improved performance in the reconstruction and prediction of instantaneous surface elevation, surface velocity potential and high-order wave statistics, including the skewness and kurtosis.


Introduction
In recent years, with the increasing capabilities in water wave measurement, substantial efforts have been made to assimilate the observation data of water waves into computational models for wave field reconstruction and prediction. In modern observational studies, the key wave properties and the spatial distribution of the wave surface elevations and velocities can be measured using remote sensing Flow E2-3 In the present study, we propose a new adjoint-based data assimilation method, named the connected-parameter method (CPM), to address the inconsistency between measurement resolution and simulation resolution. The key feature of our method is the consideration of the wave-physics-based connection between the control variables, i.e. the initial surface elevation and the initial surface velocity potential. A typical wave model imposes constraints on the time evolution of the wave state, not on the initial wave state. The unconstrained initial wave states have a significant impact on wave dynamics. Therefore, for wave reconstruction, it is necessary to develop a wave-physics-based connection to guide the algorithm to search for the optimal initial wave state. We show that the conventional method, named the free-parameter method (FPM) in this paper, has unsatisfactory performance when the measurement data have a lower resolution than the reconstructed wave field. The wave reconstruction and prediction performance of both methods are evaluated for wave data of various nonlinearity and noise levels, and the new CPM is shown to have much improved performance.

Mathematical model and methodology
In this section, we introduce the mathematical foundations of the wave reconstruction and prediction framework. As shown in figure 1, the key components include a wave model, the corresponding adjoint model and an optimiser. The wave model serves as a nonlinear function that maps the control variable, i.e. the initial wave condition, to a time series of wave simulation data. The adjoint model is used for calculating the gradients of a predefined cost function with respect to the control variables. The control variables are then updated iteratively through the optimisation process.

Wave model
For the wave simulation, we use the high-order spectral (HOS) method (Dommermuth & Yue, 1987;West et al., 1987), a phase-resolved wave model. Under the potential flow assumption, it can be shown that the wave system is uniquely determined by the quantities at the surface (Zakharov, 1968). The governing equations expanded to the third perturbation order are (Aragh & Nwogu, 2008;West et al., 1987;Yoon et al., 2015)  where = (x, y, t) is the surface elevation, = (x, y, z = , t) is the velocity potential at the water surface, z = / z(x, y, z = , t) is the surface vertical velocity, ∇ = ( / x, / y) is the gradient operator in the horizontal directions, x and y denote the horizontal coordinates, z denotes the vertical coordinate and L[ ] = −F −1 [|k|F [ ]] is a linear operator with |k| being the magnitude of the wavenumber. Here, F and F −1 denote Fourier transform and inverse Fourier transform, respectively. The boundary condition is assumed to be periodic and the spatial derivatives are calculated efficiently with the fast Fourier transform. The fourth-order Runge-Kutta method is used for time advancement of the evolution equations. The HOS method has been used extensively in wave simulations. More details on its numerical scheme and validation can be found in Mei, Stiassnie, and Yue (2018).

Adjoint model
Based on the wave model, i.e. (2.1) and (2.2), the corresponding adjoint model is (Aragh & Nwogu, 2008) (x, y, t) L-BFGS-b Figure 1. Wave field reconstruction and prediction scheme consisting of the HOS method, the adjoint model and a gradient-based optimiser.
where 1 = 1 (x, y, t) and 2 = 2 (x, y, t) are the adjoint variables. The state variables, and , which are predicted by the wave model, are stored and serve as the parameters in the adjoint model. The adjoint variables 1 and 2 are initialised as 0 at the final time instants. At each observation time instant, the difference between the predicted surface elevation obtained from the wave model and the measured data, i.e. ( − M ), is added to the adjoint variable 1 at the corresponding measured locations as 1 = 1 + ( − M ). The numerical scheme for integrating the adjoint model is the same as that for the wave model, except that the adjoint model is integrated backwards in time (i.e. with a negative time step) to obtain the adjoint variables at the initial time, which determine the gradients of the cost function with respect to the control variables, as explained in the next section.

Cost function and gradients
Searching for the optimal initial wave condition that minimises the difference between the reconstructed wave field and the measurement is a key step in wave reconstruction. In this study, we define a cost function to quantify this difference between the predicted surface elevation from the time evolution of the initial wave state and the measured surface elevation M , based on the L 2 -norm error (Aragh & Nwogu, 2008;Gronskis et al., 2013;Xu & Wei, 2016) (2.5) where N X and N Y denote the grid numbers of the measurement in the x and y coordinates, respectively, and N T denotes the number of time instants of the measurement used in the wave reconstruction process. For a given wave model, the predicted surface elevation is determined uniquely by the initial wave state used in the forward wave simulation. Therefore, J is a function of the initial conditions 0 and 0 , bounded by the wave model. The difference between our new method, CPM, and the conventional method, FPM, is summarised in table 1. In the FPM, both 0 and 0 are treated as the independent control variables to minimise

Method
Control parameters Gradient expression FPM 0 (x, y) and 0 (x, y) the cost function in the optimisation problem. In the optimisation process, 0 and 0 are updated at each iteration step using the gradient information J/ 0 and J/ 0 , respectively. In the new method, CPM, we derive a physics-based constraint that connects 0 to 0 by utilising the dispersion relation. The cost function is then determined entirely by 0 . Therefore, in an iteration step in the optimisation process of CPM, 0 is updated and 0 is calculated based on 0 via the physics-based constraint.
Here we present the expression of the gradients of the cost function with respect to the control variables. The detailed derivations are given in the supplementary material available at https://doi.org/ 10.1017/flo.2021.19. As shown in table 1, the gradients of the cost function regarding 0 and 0 in FPM The total derivative of the cost function with respect to 0 in CPM is which utilises the relation between 0 and 0 , where = (g|k|) 1/2 is the wave angular frequency and the sgn function is defined as follows: We stress that, in the CPM, the first-order approximation shown in (2.8) only applies to the initial wave field, and the nonlinear wave dynamics are captured using the nonlinear wave evolution model shown in (2.1) and (2.2).

Wave field reconstruction and prediction
With the gradient information calculated from the adjoint model as shown in table 1, the L-BFGS method (Byrd et al., 1995;Zhu et al., 1997) is then used to optimise the control parameters to reduce the cost function. Similar reconstruction frameworks, including the forward model, the adjoint model and a gradient-based optimiser, have been used in previous studies (see e.g. Aragh & Nwogu, 2008;Foures et al., 2014;Gronskis et al., 2013;Xu & Wei, 2016). The key steps to reconstruct and predict the wave field from measurement are sketched in figure 1 and summarised as follows: • Step 1. An initial guess of 0 and 0 is given for starting the wave simulation.
• Step 2. The HOS method is used for the forward simulation of the wave field from the initial time t 0 to the final time t f . • Step 3. The cost function J is calculated using (2.5). If the difference of J in two consecutive optimisation iterations is smaller than a predefined threshold value, the process ends and the initial field is the optimal solution. Otherwise, go to Step 4. • Step 4. The adjoint model is integrated from the final time t f to the initial time t 0 to obtain the gradient information J/ 0 and J/ 0 . • Step 5. The gradient information and the cost function are fed into the L-BFGS method for the optimisation of the initial condition 0 and 0 to reduce the cost function J. • Step 6: Return to Step 2. A new optimisation iteration is started with the modified initial condition 0 and 0 .

Generation of wave data
In this study, we use the wave solution obtained from the wave simulation using the third-order HOS method as the true wave data, which is sufficient to capture the nonlinear four-wave interactions (Hasselmann, 1962), to test the performance of the wave reconstruction methods. The initial condition of the wave field is constructed from the directional Joint North Sea Wave Project (JONSWAP) spectrum (Hasselmann et al., 1973) where p is the Phillips parameter, is the wave frequency, p = 1.57 rad s −1 is the peak wave frequency, p = 24.98 m is the peak wavelength, T p = 4 s is the peak wave period, = 3.3 is the peak-enhancement parameter, = 0.07 for ≤ p , = 0.09 for > p , and D( ) = 2 cos 2 ( )/π with ∈ [−π/2, π/2] is the angular spreading function. The parameters chosen here are similar to those in Qi, Wu, Liu, Kim, and Yue (2018). The computational domain size is set to L x × L y = 16 p × 16 p , with 512 grid points in the x and y directions, respectively. In each case, the wave data are collected after a relaxation period of wave evolution. In the present study, this relaxation period is 100 s, which is sufficient for capturing the nonlinear wave dynamics, as suggested in Dommermuth (2000). The simulation time interval is 0.08 s, and the time duration used for wave reconstruction and prediction is 100 s. The simulated surface elevations are referred to as the true wave field T (x, y, t).
In the simulation of wave data, we consider different wave nonlinearity and noise levels in the computational cases listed in table 2. We use two quantities to measure the wave field nonlinearity, including the effective wave steepness defined as (Qi, Wu, Liu, Kim, & Yue, 2018) is the root mean square of the initial surface elevation, and the local maximal wave steepness (ka) l = max ( 2 x + 2 y ) 1/2 . As listed in table 2, we consider a range of wave steepness in cases KA03-N00, KA06-N00, KA09-N00 and KA13-N00. To account for the effect of the measurement error in cases KA09-N03, KA09-N06 and KA09-N10, we add a random noise with a magnitude of 3 %, 6 % and 10 %, respectively, of the maximal value of wave surface elevation to the true wave field as the measurement M (x, y, t). As shown in figure 2 and table 3, M has a much lower spatial resolution than the true wave field T .
After the synthetic measurement data M are obtained, we then perform the data assimilation process separately using the FPM and CPM (see § 2.4). The grid number used for the wave reconstruction/prediction is 512 × 512. The degrees of freedom of the independent control variables is then 5.2 × 10 5 in FPM and 2.6 × 10 5 in CPM. In the optimisation iteration, we set both 0 and 0 to zero as the initial guess. The wave fields obtained by solving the governing equations (2.1) and (2.2) from the  initial guess with data assimilation using FPM and CPM are referred to as the reconstructed/predicted wave field FPM and CPM , respectively (table 3). We use the first 50 s of data for wave reconstruction and the remaining 50 s for wave prediction. Note that the wave measurement data M are assumed unknown in the prediction time duration, i.e. [50, 100] s, in the data assimilation process. The choice of the above parameters is consistent with the recent studies on wave field reconstruction and prediction (Qi, Wu, Liu, Kim, & Yue, 2018).

Performance comparison of CPM and FPM
In this section, we evaluate the performance of CPM and FPM by comparing their results with the ground truth. The results presented are from case KA09-N00, and those from other cases are presented in the next section to examine the effects of nonlinearity and measurement noise. We evaluate the convergence  Figure 3 shows the convergence of the normalised cost function as the number of optimisation iterations increases. In FPM, the cost function saturates at approximately 20 % of the initial value, while in CPM, the cost function converges to below 0.1 % of the initial value.

Wave evolution
To evaluate the algorithm performance, we first present the reconstructed initial wave state. In figure 4, we plot the instantaneous wave field at t = 0 of the true wave field T in the KA09-N00 case as described in § 3.1, the wave surface elevation FPM and CPM reconstructed using the conventional method, FPM, and our new method, CPM, respectively, and their zoomed-in views. In figure 5, we plot the same figures for velocity potential . As shown, the reconstructed surface elevation using FPM preserves the main features of the true wave field. However, we also observe spurious surface fluctuations in the wave field reconstructed by the FPM (figure 4b).
The zoomed-in view shows that the FPM overestimates the wave crests and troughs, comparing figures 4(d) and 4(e). As shown in figures 4(d) and 4( f ), the reconstructed wave field using the CPM agrees well with the true wave field, including the regions near the wave troughs and crests. The result for the surface velocity potential is shown in figure 5. The magnitude of the reconstructed velocity potential by the FPM (figure 5b,e) is significantly underestimated compared with the ground truth (figure 5a,d). On the other hand, the results produced by our CPM (figure 5c, f ) and the ground truth are indistinguishable. Next, we present the omnidirectional wavenumber spectra of reconstructed 0 of FPM, CPM and the ground truth in figure 6(a). The difference between the FPM result and the ground truth is small in the low-wavenumber range but significant in the high-wavenumber region of k > 1.5k p , while the spectrum calculated from the CPM agrees with the true value throughout the entire wavenumber range. The unsatisfactory performance of FPM is likely caused by the coarse measurement resolution, which is much smaller than the resolution of the ground-truth wave field. Specifically, for the high-wavenumber wave components that are not resolved spatially, their dynamics may still be partially captured by the time series of the wave measurement data. The FPM fails to reconstruct these wave components, while CPM is effective because of the additional constraint in (2.8). The advantage of the CPM is also seen from the spectrum of the surface velocity potential plotted in figure 6(b). In our CPM, the spectrum of the velocity potential S agrees with the ground truth.
In contrast, the FPM notably underestimates the spectrum S even for the peak wave, and the corresponding magnitude of FPM is then around half of the true value T . This might be explained by the different order of gradient magnitudes in the optimisation process (see figure S1 of the supplementary material), where J/ 0 is nearly one order of magnitude smaller than J/ 0 . In the gradient-based optimisation, the control parameters' increments are proportional to the gradient magnitudes, and thus  the cost function minimisation would place more emphasis on the parameters with large gradients, i.e. 0 , resulting in the underestimated 0 in the found suboptimal solution. However, in CPM, the control parameter 0 is connected to 0 , which would not have the same problem caused by the different order of gradient magnitudes. We also present the summed omnidirectional spectrum of the surface velocities (u, v, w)| z= at both the beginning and end of the reconstruction time duration as in figures 6(c) and 6(d), where the surface velocities are calculated from the surface elevation and surface velocity potential as (Aragh & Nwogu, 2008;West et al., 1987;Yoon et al., 2015) The spectrum in FPM deviates from the ground truth significantly and has a non-physical high energy concentration in the large-wavenumber region, whereas the result in CPM agrees well with the ground truth. Furthermore, the spectrum distribution in FPM, which changes rapidly from t = 0 to t = 50 s as evidenced by the comparison between figures 6(c) and 6(d), also indicates that FPM fails to find a physically optimal solution. We further compare the time history of the reconstructed and predicted wave fields with the true data. Considering that the reconstructed and predicted wave field has a higher spatial resolution than the measurement (see table 3), we have first examined locations where the synthetic surface elevation measurement data are available (not plotted due to space limit): the FPM generally produces a slightly worse result compared with the CPM. As a comparison, we plot in figure 7 the results obtained using both methods at a fixed location without the measurement data (x = 9.4 p , y = 8 p ). The reconstructed surface elevation obtained by FPM has notable deviations from the true wave state, and the surface velocity potential obtained by FPM shares a similar distribution with the ground truth but with a visible difference, while the results obtained by CPM agree well with the true values. In addition, we also compare the spatial distribution of the reconstructed wave field and the true wave field along the line y = 8 p at t = 24 s and t = 72 s in figure 8. As shown in figures 8(a) and 8(b), the reconstructed surface elevation using FPM contains non-physical high-wavenumber oscillations in the spatial domain. The velocity potential obtained by FPM has a noticeable difference compared to the true data T as shown in figures 8(c) and 8(d), which would result in a more significant difference in velocity as shown in figures 6(c) and 6(d) considering the fact that the velocity is related to the spatial derivative of the velocity potential. Similar to figure 7, CPM produces more accurate results compared to FPM. We also observed that, for the results computed at other locations (not plotted), the CPM always outperforms FPM. In summary, including a constraint between 0 and the independent control variable 0 in CPM provides an apparent improvement over FPM in recovering the true wave dynamics in both the reconstruction and prediction time duration.

Wave statistics
To evaluate the statistics of the reconstructed wave field, we calculate the probability density function (p.d.f.) of the wave field. In figure 9, we plot the p.d.f. of / 2 1/2 and / 2 1/2 at t = 0, t = 24 s for the reconstructed wave field, and t = 72 s for the predicted wave field, where the bracket · · · denotes the spatial mean. Under the linear wave assumptions, the p.d.f. of the surface elevation yields a Gaussian distribution. However, as observed in field measurements (see e.g. Ochi & Wang, 1985), if the wave slopes are not small, the p.d.f. deviates from the Gaussian due to the nonlinearity, which is consistent with the p.d.f. of T and T in our result. The reconstructed and predicted wave fields obtained by CPM successfully recover the p.d.f. of the true wave field, while the results by FPM have a non-negligible deviation, indicating the difference of FPM and FPM from the true values, consistent with the results in the preceding section. Skewness and kurtosis are important statistics to reflect the physical features of a nonlinear wave field. Specifically, the skewness measures the deviation of the wave profile from a sinusoidal shape, and the kurtosis indicates the probability of the occurrence of extreme waves (Xiao, Liu, Wu, & Yue, 2013). We compute the skewness C 3 and the kurtosis C 4 from the instantaneous surface elevation as (3.3a,b) We present the evolution of skewness and kurtosis of the reconstructed wave field and the true wave field in figure 10. For a standard Gaussian distribution, their values are C 3 = 0 and C 4 = 3, respectively.  x/λ p

Figure 9. Probability density function of the normalised reconstructed wave fields obtained using FPM and CPM, and the true wave field at (a,d) t = 0; (b,e) t = 24 s (in reconstruction time duration); and (c, f) t = 72 s (in prediction time duration). The standard Gaussian distribution is also plotted.
The skewness of the computed wave field varies from 0.05 to 0.2 and the kurtosis varies from 2.8 to 3.3, mostly above 3, which differs from the statistics of a standard Gaussian distribution because of the wave nonlinearity. As shown in figure 10, our CPM can successfully recover the skewness and kurtosis of the wave field in both the reconstruction time duration and the prediction time duration. On the other hand, there exists a distinct difference between the FPM results and the ground truth, especially for the kurtosis (see e.g. t/T p = 11 in figure 10b).

Effect of nonlinearity and measurement noise
To evaluate the effect of wave field nonlinearity and measurement noise, we perform the data assimilation for the reconstruction and prediction for wave fields with different nonlinearity and noise levels, as shown in table 2. The wave data generation process is described in § 3.1. We present the omnidirectional wavenumber spectra of the reconstructed initial wave field and the true wave field in figures 11 and 12 for the cases with the largest wave steepness and noise level, i.e. cases KA013-N00 and KA009-N10, respectively. Compared with the results for the case KA09-N00 (see figure 6), the deviation of the wave surface elevation obtained by CPM from the true wave data increases in the high-wavenumber range due to the high nonlinearity and noise. Nevertheless, the performance of CPM in reconstructing/predicting the surface elevation is much better than that of FPM.
To quantify the overall performance of CPM and FPM, we define the correlation coefficient between and T as   Figure 14. Same legend as in figure 13. The results for cases KA09-N00, KA09-N03,  where N X and N Y denote the grid numbers of the reconstructed and predicted wave field in the x and y coordinates, respectively. The value of ( , T ) is a measure of the accuracy of the data assimilation scheme, and a larger value corresponds to a better accuracy. When = T , ( , T ) = 1. For the surface velocity potential, the correlation coefficient is calculated using the same definition as in the above equation.
The time-averaged correlation coefficient in the reconstruction time duration (t < 50 s) and the prediction time duration (t > 50 s) for wave fields with different wave steepness and noise level are presented in figures 13 and 14. For the FPM result, the time-averaged correlation coefficients ( R , T ) and ( R , T ) in the reconstruction time duration decrease from 0.8 to 0.6 and from 0.9 to 0.7, respectively, with increasing wave steepness, and the correlation coefficients in the prediction time duration decrease from the values in the reconstruction time duration more rapidly with higher wave steepness. For the CPM result, the correlation coefficient of the reconstructed wave field with the true wave field is above 0.9 for all the cases. Figure 14 shows that the effect of noise level in the range of 0 % to 10 % has a negligible effect on the reconstruction performance. In these cases, the correlation coefficients of the reconstructed and predicted wave field obtained by CPM are notably higher than those obtained by FPM. Therefore, the advantages of CPM over FPM are unaffected by the measurement noise and higher wave nonlinearity. Simulations with different initial random wave phases are also performed for cases KA09-N00, KA09-N10 and KA13-N00. As shown in table 4, we found that the wave reconstruction and prediction performance for both FPM and CPM, quantified by correlation coefficients, varies only slightly with the initial random wave phases.

Reasons for better performance of CPM
The better performance in CPM over FPM is likely due to three factors. First, the non-physical highwavenumber wave components have been removed by using the wave-physics-based constraint for control variables in CPM. The wave model, i.e. (2.1) and (2.2), only restricts the evolution of and . However, 0 and 0 are independent variables that serve as the initial condition (Zakharov, 1968). In wave simulations, appropriate initial states are required to ensure that the wave dynamics are captured correctly. In FPM, the initial wave states are obtained via a gradient-based optimiser that does not guarantee the physical reasonableness of the found solutions for this complex system and results in a solution with non-physical high-wavenumber waves. The additional constraint in CPM, on the other hand, helps the optimisation algorithm to search for a solution with the reasonable initial wave state quantified by 0 . If the nonlinearity were to be included in the dispersion relation, the frequency of a given wave would be a function not only of the wavenumber and amplitude of itself but also of those of all other free-wave components (Wu, 2004). Therefore, it is infeasible to write a similar formula to serve as the constraint. Besides, because this linear assumption is imposed only at the initial time, and the forward wave model captures the nonlinear evolution afterwards, we do not expect a significant change in the reconstructed and predicted wave field even if a nonlinear dispersion relation can be incorporated into the constraint between 0 and 0 . Second, the additional constraint in CPM addresses the issue of the inhomogeneity in the magnitude of the gradient in the optimisation process of FPM. In the gradient-based optimisation used in FPM, the increment of the control parameters is proportional to the gradient magnitudes, and thus the optimiser tends to modify the parameters with large gradients, i.e. 0 . However, in CPM, the control parameter 0 is connected to 0 , which would not have the same problem caused by the different order of gradient magnitudes. An alternative way to examine this issue is to use the scaling strategy. Strictly speaking, the suitable value for scaling is unknown before wave reconstruction because the velocity potential is unknown and only sparse measurements of surface elevation are available. However, for the sake of the performance testing, we assume and are known and choose them as the scaling factors for and , respectively, to ensure that the normalised variables have the same order of magnitude. As shown in table 5, while the scaling strategy enhances the performance of FPM, its cost function is still 60 times larger than that in CPM, suggesting that the performance issue in FPM cannot be solved by scaling.
Third, this seemingly counter-intuitive worse performance in FPM by using an extra control variable 0 is known as the degradation problem widely observed in the optimisation of complex systems with large degrees of freedom (He, Zhang, Ren, & Sun, 2016). Specifically, when the number of independent control variables increases, the performance of the optimiser decreases counter-intuitively such that the optimiser fails to find the global optimal solution. The degradation problem is based on the observation of an abnormal increase of the cost function with adding more free control variables. Rigorous theoretical explanation is still a research topic to address in the research community. An effective method that has been adopted by the deep neural network optimisation community to alleviate the degradation problem  Figure 15. Omnidirectional wavenumber spectra of the reconstructed initial surface elevation using measurements of different temporal resolutions and the true initial surface elevation: (a) for case and (b) for a typical wave field of wave period T p = 10 s. Results are obtained using CPM.
is to add connections in different neural layers to change the system structure. Our results show that a wave-physics-based constraint introduced in the CPM can also provide an effective solution to help the optimiser when searching for the globally optimal result.

Effect of the discretisation of measurement on CPM performance
We have shown that measurement with a coarse resolution of Δx M / p = 0.5 and Δt M /T p = 0.1 is sufficient to accurately capture the high-frequency information at the present configuration using CPM. Theoretically, by assuming that the temporal resolution of the measurement is adequately high and by using the linearised wave theory, it is possible to determine the predictable zone for irregular wave fields by the maximum and minimum wave group velocities and direction spreading angle of the wave field as well as the spatial and temporal extents of the measurement Wu, 2004). However, this theoretical work is inapplicable for discretised spatial and temporal resolutions, which are typically seen in wave measurement practice. For the case KA09-N00, CPM has a fairly good performance in the reconstructed surface elevation spectrum when the temporal resolution is changed to Δt M /T p = 0.2, as shown in figure 15(a). However, the performance declines significantly when the resolution is Δt M /T p = 0.3. We also test the reconstruction performance for another wave field with p = 156 m, T p = 10 s, (ka) e = 0.08 and (ka) l = 0.33 as shown in figure 15(b). The simulation time step Δt = 0.25 s and the domain length is 16 p ×16 p with the 512×512 grid resolution. Two simulations are conducted with this wave field with measurements of different temporal discretisation Δt M /T p = 0.125 and Δt M /T p = 0.25 and the same spatial resolution Δx M / p = 0.5. Similar to the result shown above, high-wavenumber wave components are observed with Δt M /T p = 0.25, while CPM obtains good reconstruction performance for Δt M /T p = 0.125. Note that for all the cases, reconstructed wave fields obtained in CPM have higher correlation coefficients with the true wave fields than the results obtained from FPM. Therefore, the discretisation of measurement is an important factor that affects the accuracy of CPM. In real applications, we need to choose appropriate combinations of the spatial and temporal resolutions of measurement to obtain satisfactory results.

Conclusions
In this study, we have investigated the assimilation of measurement data for the wave field reconstruction and prediction based on the HOS method and its adjoint model. In our method, we quantify the difference between the reconstructed/predicted wave field and measurement as a cost function. Compared with the conventional adjoint-based wave data assimilation method, FPM, we have introduced a physical constraint on the initial wave field in the new method, CPM, which is shown to effectively reduce the cost function in both the reconstruction and prediction. In the optimisation process, we use the L-BFGSb method to minimise the cost function, and convergence can be reached within several iterations steps. The new method can be applied to situations where the wave measurement data have a low resolution.
To evaluate the performance of the new method, CPM, we have generated wave data with different wave steepness and noise levels. The gradient information calculated using the adjoint model is validated against that obtained from the finite-difference method (see the supplementary material). We have conducted numerical tests to examine the optimisation, reconstruction and prediction performance for the new CPM and the conventional FPM. In the test results, CPM shows an advantage over FPM. Using one half of the control variables in FPM, CPM proves to be more efficient in reducing the overall cost function than FPM. We have also calculated the omnidirectional wavenumber spectra of the optimal initial wave states. It is found that non-physical high-wavenumber components are generated in reconstructed surface elevation and the magnitude of the initial surface velocity potential is underestimated in the FPM results, while the CPM results are close to the ground truth, demonstrating an improved capability of wave field reconstruction and prediction. The time history and the spatial distribution of the wave states reconstructed and predicted by CPM show significantly smaller errors than FPM. In addition, the CPM can successfully predict key wave statistics, including the p.d.f., skewness and kurtosis of the wave field. We have also investigated the effect of measurement noise and wave nonlinearity, and observed a better performance of CPM over FPM in all cases.
Finally, we remark that in this study the performance of data assimilation is evaluated using the synthetic wave data obtained from simulation, which might be limited by the assumptions employed by the wave model, e.g. potential flow and periodic boundary conditions. In applications, similar tests can be conducted using wave data with realistic noise obtained from field measurements. Other effects, such as bottom topology, ambient current, wind forcing and wave breaking, are not incorporated in the present wave model but could in principle also be included with modifications (Wu, 2004). When a physical process is significant but not captured in a wave model, the error caused by the model inaccuracy would result in an inaccurate reconstructed wave field. In future studies, it would be interesting to incorporate these effects in a modified data assimilation framework for both the wave model and the adjoint model.