Data Assimilation and Parameter Identification for Water Waves Using the Nonlinear Schrödinger Equation and Physics-Informed Neural Networks

Abstract

The measurement of deep water gravity wave elevations using in situ devices, such as wave gauges, typically yields spatially sparse data due to the deployment of a limited number of costly devices. This sparsity complicates the reconstruction of the spatio-temporal extent of surface elevation and presents an ill-posed data assimilation problem, which is challenging to solve with conventional numerical techniques. To address this issue, we propose the application of a physics-informed neural network (PINN) to reconstruct physically consistent wave fields between two elevation time series measured at distinct locations within a numerical wave tank. Our method ensures this physical consistency by integrating residuals of the hydrodynamic nonlinear Schrödinger equation (NLSE) into the PINN’s loss function. We first showcase a data assimilation task by employing constant NLSE coefficients predetermined from spectral wave properties. However, due to the relatively short duration of these measurements and their possible deviation from the narrow-band assumptions inherent in the NLSE, using constant coefficients occasionally leads to poor reconstructions. To enhance this reconstruction quality, we introduce the base variables of frequency and wavenumber, from which the NLSE coefficients are determined, as additional neural network parameters that are fine tuned during PINN training. Overall, the results demonstrate the potential for real-world applications of the PINN method and represent a step toward improving the initialization of deterministic wave prediction methods.

1. Introduction

The field of ocean engineering and experimental water wave research highly desires a deterministic description of wave quantities [1], which are usually described by partial differential equations (PDEs). Unlike statistical methods, deterministic prediction involves a phase-resolved tracing of wave fields, i.e., the wave elevation $\eta$ as a function of space and time with high resolution. Unfortunately, acquiring such spatio-temporal data from experiments is often impractical: On the one hand, the reconstruction of wave surface elevations from real-world radar data is still an unresolved issue [2]. On the other hand, in situ measurement devices, such as wave gauges, measure the elevation time series at only a few selected locations due to the high operational costs. This sparsity of information leads to ill-posed inverse problems [3] when attempting to reconstruct the complete wave elevation $\eta (x,t)$ from gauge measurements ${\eta }_{\mathrm{m}}\left(t\right)$, which might be solved by numerical or machine learning methods.

Despite substantial advancements in common numerical PDE solvers, such as finite element, finite difference, and spectral methods, their applications to highly dynamic systems still incurs remarkable computational costs. These costs primarily stem from the need for fine-grained discretizations to ensure accurate solution approximations. Furthermore, conventional grid-based numerical solvers still face unresolved challenges, especially for addressing ill-posed inverse problems [4], such as wave reconstructions. In particular, inverse problems become computationally expensive due to the high number of forward evaluations required to estimate the inverse. Additionally, issues related to numerical stability and convergence necessitate the implementation of regularization techniques [5].

In the meantime, machine learning methods have revolutionized many scientific disciplines, including their applications in the fields of fast deterministic ocean wave prediction cf. [6,7,8,9] and reconstruction cf. [2,10]. However, in contrast to classical numerical solvers, these data-based approaches lack the incorporation of inherent knowledge about the physical laws that approximate the underlying system in the form of PDEs. Physical information is solely provided by observational training data for supervised learning, which are sourced from PDE simulations using classical numerical methods or from field measurements of real physical wave systems. Thus, the solutions generated by data-based approaches may not inherently ensure physical consistency, as the quality and quantity of training data limit their accuracy. This can be regarded as a neglect of established knowledge, especially when addressing the reconstruction problem outlined earlier, considering the centuries-long development of various model equations for describing surface gravity waves.

To overcome the limitations of both explicit numerical PDE solvers and neural network approaches, Raissi et al. [11] proposed physics-informed neural networks (PINNs). PINNs integrate observational data with physical laws by parameterizing the PDE solution as a neural network to solve forward, inverse parameter identification or inverse data assimilation problems. More precisely, the network’s training process is constrained by a loss function that incorporates a PDE residual. To effectively calculate this residual, the PINN algorithm leverages the method of automatic differentiation [12]. In recent years, PINNs have gained attention across diverse scientific domains [4,13], including computational fluid dynamics [14,15,16], acoustic wave propagation [17], heat transfer [18,19], climate modeling [20], nano-optics [21], and the study of hyperelastic materials [22]. Notably, PINNs also employed in a few investigations related to water waves: For instance, Wang et al. [23] used a loss function based on the wave energy balance equation and the linear dispersion relation to reconstruct near-shore phase-averaged wave heights. The same equations also allow for solving sea bed depth inversion problems from statistical wave parameters in shallow-water regimes [24]. Furthermore, the Saint-Venant equations within the loss residuals of PINNs enable the downscaling of large-scale river models by assimilating remote sensing data in conjunction with in situ measurement data of the water surface [25]. Additionally, the research of Jagtap et al. [3] showcased the potential of PINNs in resolving ill-posed assimilation problems by leveraging analytical solitary surface measurements and the Serre–Green–Naghdi equations in shallow water.

However, in typical ocean engineering research, i.e., down-scaled testing in wave tanks, scenarios arise where the water depth significantly exceeds the wavelength. In such cases, it becomes essential to characterize the nonlinear behavior of water waves in the intermediate-to-deep water regime [26] rather than in shallow water. Zakharov [27] demonstrated that the amplitude envelope of slowly modulated wave groups approximately satisfies the nonlinear Schrödinger equation (NLSE). Consequently, various variants of the hydrodynamic nonlinear Schrödinger equation were investigated experimentally and numerically for deterministic wave prediction cf. [28,29,30] and rogue wave modeling cf. [31,32].

To the best of the authors’ knowledge, there is no prior documentation of PINNs to solve the hydrodynamic form of the NLSE, despite the successful integration of other forms of the NLSE for different physical phenomena into the PINN framework cf. [11,33,34,35,36,37,38]. However, unifying most of this related work demonstrates the application of the NLSE-PINN methodology concerning the initial or boundary conditions derived from analytical soliton or breather solutions that correspond to their specific NLSE used in the loss function. While these approaches provide a theoretical proof of concept, they cannot inherently assess the practical viability of PINNs in real-world scenarios, as analytical solutions seldom align with the complexity and imperfections encountered in measurements.

Therefore, the objective of this study was to demonstrate the application of an NLSE-PINN framework in the more realistic scenario of irregular deep water gravity waves. Specifically, we aimed to leverage the capability of PINNs in solving data assimilation problems, which is crucial in ocean engineering to infer the dynamic system’s state from partial or noisy surface elevation measurements ${\eta }_{\mathrm{m}}\left(t\right)$ cf. [17,25,39,40], thereby reducing the need for costly gauges while still being able to reconstruct the fully resolved wave fields $\eta (x,t)$ in between. For this purpose, we utilized irregular and nonlinear wave group data from gauge locations several meters apart within a numerical wave tank that was simulated using the high-order spectral method (HOSM) [41].

However, as our PINN loss function was constrained by the hydrodynamic NLSE associated with limitations in accommodating arbitrary irregular sea states due to the constraints in the bandwidth, steepness, and nonlinearity [1], we investigated the utilization of a slightly misspecified PINN for assimilating wave data derived from the more realistic HOSM; the HOSM also accounts for highly nonlinear, broad-banded, or directional sea states in non-breaking scenarios [42,43]. This misspecification acknowledges the common discrepancy between mathematical models (e.g., PDEs) and measurement data in complex, real-world physical systems [44], where perfect alignment between measurements and analytical solutions is uncommon, and exact knowledge of the physical equations that describe all phenomena encountered in measurements is rare. Moreover, an ongoing challenge in applying the NLSE to irregular waves with specific bandwidths lies in accurately determining its coefficients, which are derived from ratios of the base parameters of the carrier wave frequency ${\omega }_{\mathrm{p}}$ and wavenumber $k_{\mathrm{p}}$. In particular, relatively short measurement intervals and slight deviations in the measurement data from the narrow-band assumptions of the NLSE can introduce uncertainties in the determination of ${\omega }_{\mathrm{p}}$ and $k_{\mathrm{p}}$, which affect, for example, the propagation velocities, and thus, cause errors and offsets in the surface elevation reconstructions. Therefore, our hypothesis for this study was as follows:

• Alongside the assimilation of measurement data, the NLSE-PINN will facilitate fine tuning the base parameters (${\omega }_{\mathrm{p}}$ and $k_{\mathrm{p}}$) of the NLSE coefficients. We expect this approach to enhance the reconstruction performance compared with using constant coefficients predetermined from spectral wave properties.

To address this, we first present the wave tank utilized for data generation in Section 2.1. Afterward, the hydrodynamic NLSE is introduced in Section 2.2 to subsequently develop the NLSE-PINN framework and training methodology in Section 2.3. Then, we first performed a pure data assimilation task with predetermined, constant NLSE coefficients to reconstruct wave surface envelopes from only two gauge measurements in Section 3.1. These results served as a benchmark to demonstrate the enhancement of reconstruction by incorporating the NLSE coefficient base parameters as additional tuneable PINN variables in Section 3.2. Finally, Section 4 summarizes the key findings, outlines the method limitations, and suggests potential directions for future research.

2. Method

The following subsection initiates a concise overview of the wave tank specifications and the numerical method to generate the synthetic wave elevation data. Next, the nonlinear Schrödinger equation for deep water gravity waves is introduced. Subsequently, this equation is integrated to develop a physics-informed neural network, with details on its architecture and training strategy provided in the last subsection.

2.1. Numerical Wave Tank and Measurement Data Generation

The numerical wave tank experiments were based on the dimensions of a real wave tank facility at Hamburg University of Technology. This wave tank possessed a cross-sectional area of $1.5\mathrm{m}\times 1.5\mathrm{m}$ and extended over a length of $15\mathrm{m}$ (https://www.tuhh.de/mum/en/research/facilities (accessed on 11 September 2024)) and is visualized in Figure 1a. To generate waves, a flap-type board was installed on one side of the tank, while the opposite side had a beach element to minimize wave reflections. To simulate nonlinear wave propagation inside this numerical wave tank, we employed the high-order spectral method (HOSM) of order $M=4$, as formulated by West et al. [41]. This method models highly realistic, irregular nonlinear water wave surfaces on a Cartesian coordinate system $(x,z)$, with the mean free surface located at $z=0\mathrm{m}$ with z pointing upward. Assuming a Newtonian fluid that is incompressible, inviscid, and rotational, the HOSM solves the general initial boundary-value potential flow problem that is given by the Laplace equation:

$$\begin{array}{cc}\hfill {\nabla }^2\mathrm{\Phi }={\mathrm{\Phi }}_{xx}+{\mathrm{\Phi }}_{zz}& =0\hfill \end{array}$$

(1)

and the kinematic surface, dynamic surface, and bed boundary condition:

$$\begin{array}{cc}\hfill {\eta }_t+{\eta }_x{\mathrm{\Phi }}_x-{\mathrm{\Phi }}_z& =0\mathrm{on}z=\eta (x,t)\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_t+g\eta +{\displaystyle \frac{1}{2}}\left({\mathrm{\Phi }}_x^2+{\mathrm{\Phi }}_z^2\right)& =0\mathrm{on}z=\eta (x,t)\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_z& =0\mathrm{on}z=-d.\hfill \end{array}$$

(4)

Therein, $\mathrm{\Phi }(x,z,t)$ is the velocity potential, $\eta (x,t)$ is the free surface elevation, and g is the gravity acceleration. The simulations were initialized using wave surfaces from the JONSWAP spectra [45] in the finite depth form [46].

For the data generation, we maintained a fixed water depth of $d=1\mathrm{m}$ and employed four virtual wave gauges to provide time-series measurements along the water surface at $x_{\mathrm{g}}\in \{3,4,5,6\}\mathrm{m}$, which constituted a typical depth and measurement point spacing in educational wave tanks [47,48]. The outer gauges of the numerical wave tank were designated as locations of sparse measurements ($x_{\mathrm{g},\mathrm{meas}}\in \{3,6\}\mathrm{m}$) for the PINN training, while measurements at the inner gauges ($x_{\mathrm{g},\mathrm{test}}\in \{4,5\}\mathrm{m}$) were solely used to assess the PINNs’ performance against a ground truth after the training process. We generated a temporal wave surface elevation series $\eta (x=x_{\mathrm{g}},t)$ over an interval that spanned t = 0–60$\mathrm{s}$ to ensure that the wave packets fully passed all four gauges, while no significant wave reflections had yet occurred at the channel’s end. A graphical representation of the resulting three-dimensional data structure is exemplified in Figure 1b.

As we desired wave data that covered different wave conditions and significant wave heights $H_{\mathrm{s}}$, the sea state parameters of the peak wave frequency ${\omega }_{\mathrm{p}}=2\pi /T_{\mathrm{p}}$ (with $T_{\mathrm{p}}$ being the peak period), wave steepness $ϵ=k_{\mathrm{p}}{H}_{\mathrm{s}}/2$ (with $k_{\mathrm{p}}=2\pi /L_{\mathrm{p}}$ being the peak wavenumber and $L_{\mathrm{p}}$ being the peak wavelength), and peak enhancement factor $\gamma$ were varied as follows:

$$\begin{array}{cc}\hfill {\omega }_{\mathrm{p}}& \in \{3,4,5,6,7,8,9\}\frac{\mathrm{rad}}{\mathrm{s}}\hfill \\ \hfill ϵ& \in \{0.0125,0.0250,0.0375,0.0500,0.0750,0.1000\}\hfill \\ \hfill \gamma & \in \{1,3,6\},\hfill \end{array}$$

where a higher $\gamma$ indicates a narrow-banded spectra, as shown in Figure 2. By randomly selecting initial phase shifts of the component waves, five different elevations were generated for each $ϵ$–${\omega }_{\mathrm{p}}$–$\gamma$ combination, which resulted in 630 different wave samples in total, where each individual sample had the data structure shown in Figure 1b. For details on the nonlinear wave data generation in a numerical wave tank facility using the HOSM, the reader is referred to Klein et al. [7] and Lünser et al. [43].

2.2. Hydrodynamic Nonlinear Schrödinger Equation

In addition to employing numerical methods, such as HOSM, to approximate the potential flow Equations (1)–(4), an alternative approach is the derivation of simplified solutions. By utilizing perturbation theory around the parameters of wave steepness $ϵ=k_{\mathrm{p}}a$ and relative bandwidth $\mu =\Delta k/k_{\mathrm{p}}<<1$, small-amplitude waves and a narrow bandwidth are assumed. Here, $\Delta k$ denotes the width of the wavenumber spectrum around the peak wavenumber $k_{\mathrm{p}}$. Moreover, the boundary value problem at the unknown free surface $z=\eta (x,t)$ can be approximated using a Taylor series expansion. By truncating the perturbation expansion at order $\mathcal{O}\left({ϵ}^3\right)$, the envelope equation of the nonlinear Schrödinger equation (NLSE) is derived [27,49,50] in terms of a complex wave envelope amplitude $A(x,t)=U(x,t)+iV(x,t)$ that varies slowly compared with the phase $\vartheta =k_{\mathrm{p}}x-{\omega }_{\mathrm{p}}t+\phi$ of its underlying carrier wave $\eta (x,t)$, where $\phi$ is a phase shift. The hydrodynamic NLSE in the time-like form is as follows:

$$i\left(A_x+{\displaystyle \frac{1}{c_{\mathrm{g}}}}{A}_t\right)+\delta A_{tt}+\nu {\left|A\right|}^{2}A=0$$

(5)

which finds common application in boundary value wave tank problems [51], where

$$c_{\mathrm{g}}={\displaystyle \frac{{\omega }_{\mathrm{p}}}{2k_{\mathrm{p}}}},\delta =-{\displaystyle \frac{k_{\mathrm{p}}}{{\omega }_{\mathrm{p}}^2}},\nu =-k_{\mathrm{p}}^3$$

(6)

are the NLSE coefficients. The peak frequency ${\omega }_{\mathrm{p}}$ of the carrier wave is related to the corresponding peak wavenumber $k_{\mathrm{p}}$ through the linear dispersion relation for deep water:

$${\omega }_{\mathrm{p}}=\sqrt{g{k}_{\mathrm{p}}}.$$

(7)

The first term in Equation (5) characterizes the spatial variation of the amplitude, while the second term represents the wave propagation by the group velocity $c_{\mathrm{g}}$. The third term introduces the dispersive effect, while the last term introduces the nonlinearity. In contrast to the potential flow equations, which require an approximation in the depth direction, the NLSE is regarded as the simplest nonlinear equation for deep-water wave dynamics, as it is uniquely formulated on the water surface [52]. However, while it captures the essential aspects of nonlinear water waves, the NLSE’s limitations in nonlinearity magnitude and spectral bandwidth limit its general applicability in deterministic wave prediction [1].

In practical cases, where real-valued measurement series ${\eta }_{\mathrm{m}}(x=x_{\mathrm{g}},t)$ are available, the values of ${\omega }_{\mathrm{p}}$ and $k_{\mathrm{p}}$ are commonly derived from a spectral representation via a discrete Fourier transform $F\left(\omega \right)$ of the time domain surface elevation ${\eta }_{\mathrm{m}}\left(t\right)$. We employed the advanced method of Sobey [53] and Mansard and Funke [54] by defining

$${\omega }_{\mathrm{p}}={\displaystyle \frac{\int \omega ·F{\left(\omega \right)}^{5}d\omega }{\int F{\left(\omega \right)}^{5}d\omega }},$$

(8)

which enhances the robustness compared with determining ${\omega }_{\mathrm{p}}$ as the frequency where $F\left(\omega \right)$ attains its maximum. The corresponding $k_{\mathrm{p}}$ can be derived by the dispersion relation.

The complex amplitudes for initializing the NLSE are computed as

$$A(x_{\mathrm{g}},t)=\left[\eta (x_{\mathrm{g}},t)+i\mathcal{H}\left(\eta (x_{\mathrm{g}},t)\right)\right]·exp(-i\vartheta ),$$

(9)

where $\mathcal{H}$ denotes the Hilbert transform of the measured signal [55,56] and $\vartheta =k_{\mathrm{p}}x-{\omega }_{\mathrm{p}}+\phi$ is the carrier wave’s phase. For the other way around, the equation

$$\eta (x_{\mathrm{g}},t)=\mathrm{Re}\left[A(x_{\mathrm{g}},t)·exp\left(i\vartheta \right)\right]$$

(10)

instead allows for the reverse transformation from a complex envelope back to a carrier wave elevation after the reconstruction procedure, although this transformation is not presented in this work. The relationship between the carrier wave $\eta (x_{\mathrm{g}},t)$ and complex envelope $A(x_{\mathrm{g}},t)=U(x_{\mathrm{g}},t)+iV(x_{\mathrm{g}},t)$ is illustrated in Figure 3.

2.3. Physics-Informed Neural Network for the NLSE

In light of the simplifications inherent in the hydrodynamic NLSE concerning the generalized water wave problem discussed previously, we investigated the suitability of this equation within a PINN framework for analyzing surface waves that arose from sightly different wave physics by using a phenomenon that is termed model misspecification and commonly encountered in real-world physical systems. Our selection of this equation was motivated by the NLSE’s simple definition (Equation (5)), which is limited to the $(x,t)$ domain only. In contrast, for example, employing a PINN to solve the fully nonlinear potential flow equations (Equations (1)–(4)) would demand substantially greater human and computational effort. Hence, subsequent sections introduce the architecture, loss function, training methodology, and evaluation metrics adopted for our NLSE-PINN framework.

2.3.1. PINN Architecture

The neural network architecture of the PINN developed in our study is depicted on the lower-left side of Figure 4. It includes two input nodes for collocation points in space x and time t. Given the complex-valued solution $A(x,t)=U(x,t)+iV(x,t)$ of the NLSE, two output nodes are required [11], where U represents the real part and V represents the imaginary part. Moreover, the network consists of four intermediate hidden layers, with each comprising 200 nodes, which results in a total network depth of $D=6$. The variables of each layer l, which are denoted as ${\mathit{\theta }}^{\left(l\right)}={\{{\mathbf{W}}^{\left(l\right)},{\mathbf{b}}^{\left(l\right)}\}}_{1\le l\le D}$, receive their initial values using the Xavier initialization technique for weights ${\mathbf{W}}^{\left(l\right)}$ [57], while the biases ${\mathbf{b}}^{\left(l\right)}$ are initialized to zero. The total amount of the network’s weights and biases is denoted as $\mathit{\theta }=\left[{\mathit{\theta }}^{\left(1\right)},\dots ,{\mathit{\theta }}^{\left(D\right)}\right]$.

To enhance the convergence and accuracy during the training, we incorporated a strategy using layer-wise locally adaptive activation functions [58,59], which has shown promising applications, e.g., in Pu et al. [35], Jagtap et al. [3,60], Shukla et al. [61], and Guo et al. [62]. The input of a hidden layer ${\mathbf{h}}^{(l+1)}$ is derived from the output of the preceding layer ${\mathbf{o}}^{\left(l\right)}$ by following the rule

$$\begin{array}{cc}\hfill {\mathbf{o}}^{\left(l\right)}& ={\mathbf{W}}^{\left(l\right)\mathrm{T}}·{\mathbf{h}}^{\left(l\right)}+{\mathbf{b}}^{\left(l\right)}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\mathbf{h}}^{(l+1)}& =tanh(s·a^{\left(l\right)}·{\mathbf{o}}^{\left(l\right)})\hfill \end{array}$$

(12)

for $l=1,\dots ,D-1$. The last layer D has a linear activation function. In this context, ${\mathbf{W}}^{\left(l\right)}$ and ${\mathbf{b}}^{\left(l\right)}$ are the weight and bias matrices of the $\left(l\right)$-th layer, and $s=10$ is a fixed scaling factor. Furthermore, $\mathit{a}=\left[a^{\left(1\right)},\dots ,a^{(D-1)}\right]$ represent additional variables of the network, which are fine tuned during the training process to modulate the slopes of the activations. In our case, we initialized the system with $a^{\left(l\right)}=0.2$, which resulted in the initial slope of $s·a^{\left(l\right)}=2$ being slightly steeper than the usual tanh activation. Empirical evidence suggested that this initialization could be considered valid according to [58], as it accelerated the convergence rate in our setup without causing divergence or increasing oscillations in the loss function.

Moreover, over the course of our study (Section 3.2), we examined how the data-driven fine tuning of the NLSE coefficients (Equation (6)) while assimilating sparse measurement data affected the reconstruction quality. Given that all the coefficients in the NLSE were related to the peak wavenumber and frequency, we selected treating the coefficients’ base parameters ${\omega }_{\mathrm{p}}$ and $k_{\mathrm{p}}$ as additional PINN variables while sharing an optimizer with the neural network parameters $\mathit{\theta }$ and $\mathit{a}$. This is visualized in the lower-left box of Figure 4.

2.3.2. PINN Loss Function and Training

To enable the PINN’s capability to approximate an NLSE surrogate model $\tilde{A}(x,t)=\tilde{U}(x,t)+i\tilde{V}$, we used wave elevation measurement data ${\eta }_{\mathrm{m}}\in {\mathbb{R}}^{N_{\mathrm{d}}}$ from the two outer gauge positions, $x_{\mathrm{g},\mathrm{meas}}\in \{3,6\}\mathrm{m}$, within the numerical wave tank shown in Figure 1. The elevation data at the remaining locations, $x_{\mathrm{g},\mathrm{test}}\in \{4,5\}\mathrm{m}$, were reserved for later evaluation and not incorporated in the PINN’s training process. Each elevation measurement was transferred to a complex envelope $A_{\mathrm{m}}\in {\mathbb{C}}^{N_{\mathrm{d}}}$ using Equation (9), as depicted in the top boxes of Figure 4. In this work, $N_{\mathrm{d}}=1200$ was the total number of measurement data points at the domain boundaries, which is denoted as ${\{x_{\mathrm{d}}=x_{\mathrm{g},\mathrm{meas}},t_{\mathrm{d}}^{\left(j\right)}\}}_{j=1}^{N_{\mathrm{d}}}$. These points were obtained by sampling the temporal sequence t = 0–60$\mathrm{s}$ with an increment of $\Delta t=0.05\mathrm{s}$. In addition, we incorporated a set of $N_{\mathrm{r}}=20,000$ randomly located collocation points ${\{x_{\mathrm{r}}^{\left(j\right)},t_{\mathrm{r}}^{\left(j\right)}\}}_{j=1}^{N_{\mathrm{r}}}$ to enforce the NLSE solution across the entire computational $(x,t)$ domain. Using these measurements and sets of points, the multi-objective PINN loss function was defined by

$$\begin{array}{cc}\hfill \mathcal{L}& ={\mathcal{L}}_{U,\mathrm{data}}+{\mathcal{L}}_{V,\mathrm{data}}+{\mathcal{L}}_{U,\mathrm{res}}+{\mathcal{L}}_{V,\mathrm{res}},\hfill \end{array}$$

(13)

where the loss components were composed of the mean-squared errors (MSEs):

$$\begin{array}{cc}\hfill {\mathcal{L}}_{U,\mathrm{data}}& =\underset{{\mathrm{MSE}}_{U,\mathrm{data}}}{\underbrace{{\displaystyle \frac{1}{N_{\mathrm{d}}}}\sum _{j=0}^{N_{\mathrm{d}}}{\left|\tilde{U}(x_{\mathrm{d}},t_{\mathrm{d}}^{\left(j\right)})-U_{\mathrm{m}}^{\left(j\right)}\right|}^2}}·{\left({\lambda }_{\mathrm{d}}^{\left(j\right)}\right)}^2\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{V,\mathrm{data}}& =\underset{{\mathrm{MSE}}_{V,\mathrm{data}}}{\underbrace{{\displaystyle \frac{1}{N_{\mathrm{d}}}}\sum _{j=0}^{N_{\mathrm{d}}}{\left|\tilde{V}(x_{\mathrm{d}},t_{\mathrm{d}}^{\left(j\right)})-V_{\mathrm{m}}^{\left(j\right)}\right|}^2}}·{\left({\mu }_{\mathrm{d}}^{\left(j\right)}\right)}^2\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{U,\mathrm{res}}& =\underset{{\mathrm{MSE}}_{U,\mathrm{res}}}{\underbrace{{\displaystyle \frac{1}{N_{\mathrm{r}}}}\sum _{j=0}^{N_{\mathrm{r}}}{\left|{\tilde{R}}_U(x_{\mathrm{r}}^{\left(j\right)},t_{\mathrm{r}}^{\left(j\right)})\right|}^2}}·{\left({\lambda }_{\mathrm{r}}^{\left(j\right)}\right)}^2\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\mathcal{L}}_{V,\mathrm{res}}& =\underset{{\mathrm{MSE}}_{V,\mathrm{res}}}{\underbrace{{\displaystyle \frac{1}{N_{\mathrm{r}}}}\sum _{j=0}^{N_{\mathrm{r}}}{\left|{\tilde{R}}_V(x_{\mathrm{r}}^{\left(j\right)},t_{\mathrm{r}}^{\left(j\right)})\right|}^2}}·{\left({\mu }_{\mathrm{r}}^{\left(j\right)}\right)}^2.\hfill \end{array}$$

(17)

The terms ${\mathrm{MSE}}_{U,\mathrm{data}}$ and ${\mathrm{MSE}}_{V,\mathrm{data}}$ served to quantify the error between PINN predictions $(\tilde{U},\tilde{V})$ and measured envelope data $(U_{\mathrm{m}},V_{\mathrm{m}})$ at the domain boundaries ${x,\mathrm{meas}}_{\mathrm{g}}\in \{3,6\}\mathrm{m}$. In contrast, ${\mathrm{MSE}}_{U,\mathrm{res}}$ and ${\mathrm{MSE}}_{V,\mathrm{res}}$ quantified the degree to which the PINN solution conformed to the residuals of the NLSE, which are given by

$$\begin{array}{cc}\hfill {\tilde{R}}_U(x,t):=& -{\tilde{V}}_x(x,t)-{\displaystyle \frac{1}{c_g}}{\tilde{V}}_t(x,t)+\delta {\tilde{U}}_{tt}(x,t)+\nu \left(\tilde{U}{(x,t)}^2+\tilde{V}{(x,t)}^2\right)\tilde{U}(x,t)\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\tilde{R}}_V(x,t):=& {\tilde{U}}_x(x,t)+{\displaystyle \frac{1}{c_g}}{\tilde{U}}_t(x,t)+\delta {\tilde{V}}_{tt}(x,t)+\nu \left(\tilde{U}{(x,t)}^2+\tilde{V}{(x,t)}^2\right)\tilde{V}(x,t)\hfill \end{array}$$

(19)

within the remaining computational domain. Thus, as the optimizer minimizes these residuals during training, the PINN solution increasingly aligns with the physical constraints imposed by the NLSE. The derivative terms in these residuals are calculated using automatic differentiation (AD) [12] of the neural network’s outputs $\tilde{U}$ and $\tilde{V}$ with respect to the input variables x or t, as illustrated on the right of Figure 4. Unlike numerical differentiation methods, AD provides exact derivatives without any approximation error [13]. Furthermore, AD enhances the computational efficiency by decomposing complex functions into elementary operations and reusing intermediate results through the application of the chain rule, thus avoiding redundant calculations.

Furthermore, following the approach proposed by McClenny and Braga-Neto [63], trainable self-adaptation weights ${\lambda }_{\mathrm{d}}^{\left(j\right)}$, ${\mu }_{\mathrm{d}}^{\left(j\right)}$, ${\lambda }_{\mathrm{r}}^{\left(j\right)}$, and ${\mu }_{\mathrm{r}}^{\left(j\right)}$ were introduced in the loss components (Equations (14)–(17)). These variables, which are initially set to one, are associated with each measurement point ${\{x_{\mathrm{d}},t_{\mathrm{d}}^{\left(j\right)}\}}_{j=1}^{N_{\mathrm{d}}}$ or collocation point ${\{x_{\mathrm{r}}^{\left(j\right)},t_{\mathrm{r}}^{\left(j\right)}\}}_{j=1}^{N_{\mathrm{r}}}$. The PINN autonomously identifies challenging regions inside the solution domain characterized by high point-specific errors and increases the respective self-adaptation weights to emphasize the penalty, and thus, improve the approximation. This behavior is achieved through the concurrent minimization of the loss function $\mathcal{L}=\mathcal{L}(\mathit{\theta },\mathit{a},{\mathit{\lambda }}_{\mathrm{d}},{\mathit{\mu }}_{\mathrm{d}},{\mathit{\lambda }}_{\mathrm{r}},{\mathit{\mu }}_{\mathrm{r}})$ by updating the network weights and biases $\mathit{\theta }$, activation slopes $\mathit{a}$, and coefficient parameters ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ alongside the maximization of the loss for the self-adaption weights ${\mathit{\lambda }}_{\mathrm{d}}$, ${\mathit{\mu }}_{\mathrm{d}}$, ${\mathit{\lambda }}_{\mathrm{r}}$, and ${\mathit{\mu }}_{\mathrm{r}}$ in each epoch of training.

In many instances, PINNs are trained using a two-step strategy that involves the Adam optimizer [64] for a defined number of epochs, followed by the L-BFGS optimizer [65]. This strategy is recognized as pivotal [66], with Adam initially preventing convergence to local minima, while L-BFGS refines small-scale solution components [13]. However, this approach, which is beneficial for smooth analytical solutions, proved less suitable for our studies: The L-BFGS optimizer tends to over-refine noisy elements within measurements. We found that the AMSGrad modification of the Adam optimizer [67] from the PyTorch library (Version 1.8.1, sourced from Hamburg, Germany) [68] yielded satisfactory results with a learning rate of $\alpha =0.0001$ over $15,000$ epochs of training.

2.3.3. Evaluation

After assimilating the measurements at the domain boundaries to reconstruct the envelopes in the entire computational domain, the evaluation of the PINNs’ performance necessitated a metric that compared the measurement time series $A_{\mathrm{m}}=A_{\mathrm{m}}(x=x_{\mathrm{g}},t)$ with the reconstruction $\tilde{A}=\tilde{A}(x=x_{\mathrm{g}},t)$ at all gauge positions $x_{\mathrm{g}}\in \{3,4,5,6\}\mathrm{m}$. While metrics based on Euclidean distances, such as the mean-squared error (MSE), are scale dependent and treat deviations in frequency and phase as amplitude errors [69], the surface similarity parameter (SSP)

$$\mathrm{SSP}(A_{\mathrm{m}},\tilde{A})={\displaystyle \frac{\sqrt{\int |F_{A_{\mathrm{m}}}\left(\omega \right)-F_{\tilde{A}}{\left(\omega \right)|}^{2}d\omega }}{\sqrt{\int |F_{A_{\mathrm{m}}}{\left(\omega \right)|}^{2}d\omega }+\sqrt{\int |F_{\tilde{A}}{\left(\omega \right)|}^{2}d\omega }}}\in [0,1]$$

(20)

proposed by Perlin and Bustamante [70] combines phase, amplitude, and frequency errors into a scalar unified measure. In this metric, $\omega$ denotes the wave frequency vector and $F_{A_{\mathrm{m}}}$ and $F_{\tilde{A}}$ denote the discrete Fourier transforms of the time series $A_{\mathrm{m}}(x=x_{\mathrm{g}},t)$ or $\tilde{A}(x=x_{\mathrm{g}},t)$. The SSP is a normalized error metric, with $\mathrm{SSP}=0$ indicating perfect agreement and $\mathrm{SSP}=1$ implying a comparison against zero or phase-inverted surfaces. To visually illustrate the impact of the SSP metric, we systematically compared an exemplary wave envelope against various signals across a range of SSP and MSE values in the Appendix (Figure A1). Due to its straightforward error assessment and the applicability for comparing signals with differing sampling rates and lengths, the SSP has found application in recent research related to ocean wave prediction and reconstruction by Klein et al. [1,7], Wedler et al. [9,69], Desmars et al. [71,72], Lünser et al. [43], Kim et al. [73], and Ehlers et al. [2].

3. Results and Discussion

After a successful validation of the NLSE-PINN against an analytic solution to the NLSE (see Appendix A.2), we assessed its applicability to more realistic scenarios that involved data that do not entirely conform to the equation in the loss function, such as elevation measurements from the wave tank shown in Figure 1. The NLSE-PINN framework underwent training for all 630 generated wave measurement samples, where each was characterized by a distinct wave parameter combination and random phases. Our aim was to reconstruct wave envelopes throughout the spatio-temporal domain $3\mathrm{m}\le x\le 6\mathrm{m}$ and $0\mathrm{s}\le t\le 60\mathrm{s}$ solely based on two time series measured at the gauge locations ${x,\mathrm{meas}}_{\mathrm{g}}\in \{3,6\}\mathrm{m}$ spaced 3 m apart. Additional measurements at locations ${x,\mathrm{test}}_{\mathrm{g}}\in \{4,5\}\mathrm{m}$ were reserved exclusively for evaluating the PINN solutions post-training. Therefore, the first subsection (Section 3.1) focuses on the pure data assimilation task, where the NLSE coefficients in the PINN’s loss function remained constant during training and determined a priori based on spectral wave properties and linear dispersion. In contrast, the second subsection (Section 3.2) explores treating the peak frequency and wavenumber, which form the NLSE coefficients, as identifiable variables that are adapted by the optimizer to fit the measurements best.

3.1. Data Assimilation with Constant NLSE Coefficients

In the following benchmark data assimilation task, the NLSE coefficients (Equation (6)) were kept constant during the PINN training and determined a priori based on the mean of the spectral peak frequencies ${\omega }_{\mathrm{p},3}$ and ${\omega }_{\mathrm{p},6}$ at the domain boundaries using Equation (8). This peak frequency ${\overline{\omega }}_{\mathrm{p}}=\frac{1}{2}\left({\omega }_{\mathrm{p},3}+{\omega }_{\mathrm{p},6}\right)$ was calculated for each of the 630 samples, where each individual sample was generated following the procedure outlined in Section 2.1 and followed the data structure shown in Figure 1b. The linear dispersion relation Equation (7) yielded the corresponding peak wavenumber ${\overline{k}}_{\mathrm{p}}$. Note that the calculated ${\overline{\omega }}_{\mathrm{p}}$ could slightly deviate from the peak frequency ${\omega }_{\mathrm{p}}$ used to generate the JONSWAP spectrum for the HOSM. Each sample underwent individual PINN training for 15,000 epochs, which required approximately $1200\mathrm{s}$ of computational time on a NVIDIA GeForce RTX 3090 GPU.

Figure 5 depicts a representative training loss curve for the NLSE-PINN framework. During the initial epochs, a remarkable decrease in the PDE residual error components ${\mathrm{MSE}}_{U,\mathrm{res}}$ and ${\mathrm{MSE}}_{V,\mathrm{res}}$ was observed, while the data errors ${\mathrm{MSE}}_{U,\mathrm{data}}$ and ${\mathrm{MSE}}_{V,\mathrm{data}}$ remained high. As the data errors started to decrease and the PINN solution aligned more closely with the prescribed boundary data, the PDE errors experienced a transient increase before reaching a plateau. These observations could have arose from a zero envelope that caused high data errors in the early training, despite satisfying PDE residues that caused low PDE errors. However, as the envelope progressively matched the actual measurement data at the boundaries, fulfilling PDE residues in the remaining domain became more challenging, which led to a slight increase in the corresponding error.

Figure 6 presents a reconstruction of the real part $\tilde{U}(x,t)$ and imaginary part $\tilde{V}(x,t)$ of an envelope above the carrier wave sample no. 1, with $ϵ=0.075$, ${\overline{\omega }}_{\mathrm{p}}=4.45{\displaystyle \frac{\mathrm{rad}}{\mathrm{s}}}$, and $\gamma =6$. The PINN reconstructed a wave envelope structure across the entire computational domain, despite relying solely on measurement data from the boundaries marked in violet.

While some wave samples were successfully assimilated, challenges were evident for other samples. For instance, Figure 7 presents a PINN solution for sample no. 2, which was characterized by $ϵ=0.025$ and ${\overline{\omega }}_{\mathrm{p}}=8.34{\displaystyle \frac{\mathrm{rad}}{\mathrm{s}}}$, with a broad-banded nature indicated by $\gamma =1$. The average errors for the real and imaginary parts of the envelope were $\mathrm{SSP}=0.194$ and $\mathrm{SSP}=0.242$, respectively. Notably, increased errors were observed within the computational domain at ${x,\mathrm{test}}_{\mathrm{g}}\in \{4,5\}\mathrm{m}$, despite the satisfactory alignment of the boundary values with the ground truth at ${x,\mathrm{meas}}_{\mathrm{g}}\in \{3,6\}\mathrm{m}$. These errors seemed to originate from nonphysical envelope peaks that emerged from both boundaries toward the middle of the domain but met with an offset, which indicated a small error in the NLSE coefficient of the group velocity $c_{\mathrm{g}}=\frac{{\overline{\omega }}_{\mathrm{p}}}{2{\overline{k}}_{\mathrm{p}}}$ for the underlying data. However, as the advanced method of Equation (8) was already utilized for ${\overline{\omega }}_{\mathrm{p}}$ determination, it became apparent that establishing a robust approach to determine this parameter consistently across all samples presented a challenge.

The trends observed in the aforementioned examples were repetitive across all 630 instances in the assimilation task, which is summarized in Figure 8, where the SSP value represents the average error across the real and imaginary parts of each sample. Each cell in the figure corresponds to the mean SSP achieved for all five samples of a specific $ϵ$–${\overline{\omega }}_{\mathrm{p}}$–$\gamma$ combination, where the determined ${\overline{\omega }}_{\mathrm{p}}$ value was rounded to the nearest integer for this illustration. Notably, the errors tended to decrease as the ${\overline{\omega }}_{\mathrm{p}}$ value increased. The frequency dependence aligned with the inherent preference of PINNs for successive solving from lower to higher frequency components cf. [13] when considering that higher frequencies ${\omega }_{\mathrm{p}}$ of the carrier wave $\eta (x,t)$ resulted in smoother, lower-frequency envelopes $A(x,t)$, while lower frequencies ${\omega }_{\mathrm{p}}$ of the carrier wave $\eta (x,t)$ resulted in higher-frequency envelopes $A(x,t)$. Moreover, Figure 8 illustrates a general error increase as the peak enhancement factors $\gamma$ decrease or the steepness values $ϵ$ slightly increased. As proven by Appendix A.2, these observations were not attributable to the PINN method itself but rather aligned with the model misspecification. Specifically, the NLSE was limited for narrow-band and small-amplitude waves, while the ground truth HOSM data seem to exceed its validity range for increasing $ϵ$ and a broader spectra.

3.2. Coefficient Fine Tuning Alongside Data Assimilation

Despite employing an advanced formula (Equation (8)) to determine the peak frequency ${\overline{\omega }}_{\mathrm{p}}$ required for calculating the NLSE coefficients (Equation (6)), uncertainties in the ${\overline{\omega }}_{\mathrm{p}}$ determination persisted due to the short measurement durations and broader-banded signals, as discussed in the context of the envelope offset observed in Figure 7. Given the impracticality of exploring alternative possibilities for ${\overline{\omega }}_{\mathrm{p}}$ determination or iteratively adjusting it to recalculate the solution, we examined a data-driven opportunity to fine tune this parameter in the following subsection. Therefore, we treated the PDE coefficients not as constants but as additional trainable variables to address the inverse problems during the assimilation of the measurement data. Our objective was to enhance the reconstruction quality and eliminate the envelope offsets attributed to non-optimal constant coefficients by utilizing the results from the previous Section 3.1 as a benchmark. Since all the NLSE coefficients depended on ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$, we refrained from directly setting $c_{\mathrm{g}}$, $\delta$, and $\nu$ as variables. Instead, we initially determined ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ using Equations (7) and (8), as demonstrated before, but subsequently treated them as additional PINN variables. Hence, the coefficients $c_{\mathrm{g}}$, $\delta$, and $\nu$ emerged from the current ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ in each epoch of training. Moreover, as water waves with increasing nonlinearity are known to deviate from the linear dispersion relation (Equation (7)), we decoupled the current ${\overline{k}}_{\mathrm{p}}$ and ${\overline{\omega }}_{\mathrm{p}}$ during training. As discussed in the previous subsection, each of the 630 samples, with a data structure shown in Figure 1b, underwent individual training for 15,000 epochs. However, now the training incorporated two additional trainable variables ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$, as illustrated in the lower-left box of Figure 4.

The visual comparison of the new reconstruction for sample no. 2 that employed trainable variables ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ in Figure 9 with the previous solution that used constant coefficients in Figure 7 revealed a strong reduction in the envelope offset.

Table 1. NLSE coefficients of sample no. 2 calculated from ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ for the pure data assimilation task (initial) compared with the optimized values learned during the parameter fine tuning. In addition, the corresponding SSP and MSE errors are given.

			NLSE Coefficients			Mean Error
	${\overline{\mathit{\omega }}}_{\mathbf{p}}$	${\overline{\mathit{k}}}_{\mathbf{p}}$	${\mathit{c}}_{\mathbf{g}}=\frac{{\overline{\mathit{\omega }}}_{\mathbf{p}}}{\mathbf{2}{\overline{\mathit{k}}}_{\mathbf{p}}}$	$\mathit{\delta }=-\frac{{\overline{\mathit{k}}}_{\mathbf{p}}}{{\overline{\mathit{\omega }}}_{\mathbf{p}}^{\mathbf{2}}}$	$\mathit{\nu }=-{\overline{\mathit{k}}}_{\mathbf{p}}^{\mathbf{3}}$	$\mathbf{SSP}$  Re.	$\mathbf{SSP}$  Im.	$\mathbf{MSE}$  Re.	$\mathbf{MSE}$  Im.
Initial	8.344	7.097	0.588	−0.102	−357.4	0.194	0.242	$2.72$ × ${10}^{-7}$	$2.88$ × ${10}^{-7}$
Learned	9.355	6.202	0.754	−0.071	−256.7	0.150	0.188	$1.72$ × ${10}^{-7}$	$1.73$ × ${10}^{-7}$

compares the initial peak frequency ${\overline{\omega }}_{\mathrm{p}}$ and wavenumber ${\overline{k}}_{\mathrm{p}}$ with their fine tuned values, along with the resulting NLSE coefficients. For this sample, the ${\overline{\omega }}_{\mathrm{p}}$ increased during the training, while the ${\overline{k}}_{\mathrm{p}}$ decreased, which resulted in a faster group velocity $c_{\mathrm{g}}=\frac{{\overline{\omega }}_{\mathrm{p}}}{2{\overline{k}}_{\mathrm{p}}}$. Simultaneously, the nonlinear term received slightly less attention due to a decreased coefficient $\nu =-{\overline{k}}_{\mathrm{p}}^3$. These adjustments not only resulted in a visual reduction of the nonphysical envelope offsets but also improved the SSP error values. For the real part, the error decreased from $\mathrm{SSP}=0.194$ to $\mathrm{SSP}=0.150$, and for the imaginary part, from $\mathrm{SSP}=0.242$ to $\mathrm{SSP}=0.188$, both of which were an around $20\%$ improvement for sample no. 2. This improvement is noteworthy considering that while adjustable coefficients can mitigate uncertainties in ${\overline{\omega }}_{\mathrm{p}}$ determination, they cannot fully overcome the NLSE’s assumption of a single ideal carrier frequency, particularly for HOSM data with $\gamma =1$.

Examining the 630 samples collectively revealed that the solutions that used the trainable ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ values exhibited fewer envelope offsets. Figure 10 demonstrates an enhancement in the SSP values when comparing these results for fine tuned coefficients with results obtained with the constant and predetermined coefficients shown in Figure 8. While the trend persisted, where lower-frequency, higher-steepness, or broad-banded samples yielded comparatively high individual errors, the average SSP values experienced a reduction from $\mathrm{SSP}=0.223$ to $\mathrm{SSP}=0.209$ ($\gamma =1$), from $\mathrm{SSP}=0.199$ to $\mathrm{SSP}=0.185$ ($\gamma =3$), and from $\mathrm{SSP}=0.181$ to $\mathrm{SSP}=0.165$ ($\gamma =6$). On average the reduction from $\mathrm{SSP}=0.204$ with constant coefficients to $\mathrm{SSP}=0.186$ with trainable coefficients represented an improvement of approximately $8.8\%$, while complicated instances, such as sample no. 2, also experienced considerably higher individual improvements (around $20\%$). Consequently, allowing the fine tuning of ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ values during PINN training confirmed the initial hypothsis, as it enhanced the measurable and visible reconstruction quality compared with using constant coefficients in pure data assimilation tasks.

During the coefficient fine tuning, the peak frequency increased and the wavenumber decreased while training of sample no. 2 shown in Figure 9 and Table 1. In contrast, other samples displayed varied trends, such as concurrent increasing or decreasing ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$, or decreasing ${\overline{\omega }}_{\mathrm{p}}$ with increasing ${\overline{k}}_{\mathrm{p}}$. Figure 11 illustrates this variability across all 630 samples according to their respective $ϵ$–$\gamma$ combination. The initial ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ values aligned with the blue graph of the linear dispersion relation, while the fine tuned values mostly deviated from this relation, as indicated by red crosses. The left and middle columns of this Figure reveal that medium- and broad-banded sea states ($\gamma =3$ and $\gamma =1$) tended to benefit from ${\overline{\omega }}_{\mathrm{p}}$–${\overline{k}}_{\mathrm{p}}$ combinations above the linear dispersion relation. This suggests that increased group velocities $c_{\mathrm{g}}=\frac{{\overline{\omega }}_{\mathrm{p}}}{2{\overline{k}}_{\mathrm{p}}}$ enhanced the reconstruction quality. In contrast, narrow-banded samples ($\gamma =6$) showed learned combinations both above and below the linear dispersion curve, which reflected that the neural network effectively reduced the uncertainties in the initial determination of ${\overline{\omega }}_{\mathrm{p}}$.

Figure 11 further illustrates a consistent trend: irrespective of the $\gamma$ value, increasing steepness $ϵ$ resulted in smaller changes in the ${\overline{k}}_{\mathrm{p}}$ values during training. This was evident from the vertically oriented lines between the initial and final ${\overline{\omega }}_{\mathrm{p}}$–${\overline{k}}_{\mathrm{p}}$ combinations in the upper subplots. This observation implies that as the steepness increases, the PINN prioritizes maintaining the nonlinear term of the NLSE, as represented by the coefficient $\nu =-{\overline{k}}_{\mathrm{p}}^3$. This aligns with physical reasoning, as a higher $ϵ$ corresponds to increased nonlinearities. Preserving this coefficient indicates that the PINN identifies nonlinearity as crucial. Consequently, alterations in the group velocity $c_{\mathrm{g}}=\frac{{\overline{\omega }}_{\mathrm{p}}}{2{\overline{k}}_{\mathrm{p}}}$ primarily arise from alterations in ${\overline{\omega }}_{\mathrm{p}}$ only.

4. Conclusions

This study demonstrated realistic application scenarios of physics-informed neural networks (PINNs) in ocean engineering contexts. Specifically, the hydrodynamic nonlinear Schrödinger equation (NLSE) constrains the developed PINN framework to ensure the physical consistency of its solutions. First, the NLSE-PINN is employed for reconstructing synthetic wave elevation measurements between two gauges, a process known as data assimilation and a fundamental function of PINNs. However, the utilized wave data occasionally deviates from the small-steepness and narrow-band assumptions inherent to the NLSE, which poses a challenge known as model misspecification. This challenge must be anticipated whenever employing PINNs for real-world scenarios, where analytical solutions to the PDE utilized for a PINN’s loss functions are typically not observable, and comprehensive knowledge of physical equations that cover all observed phenomena is rare. Moreover, this model misspecification and the short time intervals of elevation measurements complicate the determination of the required NLSE coefficients based on spectral properties. These uncertainties in the a priori determination of constant NLSE coefficients occasionally lead to non-physical envelope offsets in the reconstruction of certain instances.

To enhance the reconstruction quality of such cases, conventional numerical methods would typically require a systematic coefficient adjustment and numerical forward propagation to optimize the data fit between two measurement points. In contrast, our second research phase demonstrated that the NLSE-PINN method allowed for fine tuning the NLSE coefficients’ base parameters ${\overline{\omega }}_{\mathrm{p}}$ and ${\overline{k}}_{\mathrm{p}}$ concurrently with the data assimilation in the same optimization process. This integrated parameter fine tuning effectively mitigated the previously observed non-physical offsets and reduced the mean error across all 630 wave samples from $\mathrm{SSP}=0.204$ with predetermined, constant coefficients during pure data assimilation to $\mathrm{SSP}=0.186$ with trainable coefficients. This represented an average improvement of $8.82\%$, while some samples experienced even greater improvements of up to $20\%$.

However, it is important to note that the reconstruction remains generally more challenging for broader wave spectra and samples with higher steepness. This challenge is not inherent to the PINN methodology itself, but rather associated with the characteristics of the utilized NLSE to constrain the loss function for data stemming from the more realistic HOSM. Therefore, future research can explore the integration of the modified nonlinear Schrödinger equation (MNLSE), which is also referred to as the Dysthe equation [74], into the PINN loss function. Due to its higher order of nonlinearity, the MNLSE may mitigate the steepness limitations inherent in the standard NLSE and provide new insights. Also, the broader-bandwidth variant of the MNLSE [75] might be particularly advantageous, as it additionally alleviates the limitations associated with narrow-band assumptions. Moreover, the development of a PINN capable of directly solving the fully nonlinear potential flow equations represents a promising direction, as it could offer improvements over relying on simplified models, such as the variants of the NLSE.

Moreover, concerning the successful fine tuning of $\overline{{\omega }_{\mathrm{p}}}$ and ${\overline{k}}_{\mathrm{p}}$ values in this study, further research can explore the feasibility of deriving a nonlinear dispersion relation as a function of wave parameters, such as the bandwidth, amplitude, or steepness. This research direction holds the potential to deepen our understanding of the interplay of NLSE coefficients and wave parameters and may contribute to the development of more flexible and accurate models for describing nonlinear water wave phenomena. Furthermore, this study employed synthetic wave data generated at two locations within a numerical wave tank with a fixed measurement distance. Hence, future research should investigate the impact of varying the distances between these two buoys. Moreover, applying the potential flow PINN or NLSE-PINN methods to real wave measurement data from wave tanks or open ocean scenarios, particularly for directional two-dimensional wave surfaces, would be valuable.