A machine learning architecture for including wave breaking in envelope-type wave models

.


Introduction
Surface gravity wave breaking is a familiar phenomenon to those observing the ocean (Babanin, 2011), yet a good understanding capturing the fundamental physics of wave breaking has proved difficult for scientists and engineers.Understanding breaking is important for multiple reasons, including for understanding the energy balance in the ocean (Hasselmann, 1974), ocean mixing (Melville et al., 1998), airsea interaction (Melville, 1996;Deike, 2022), short-term wave statistics (Toffoli et al., 2010), and loading on maritime structures (Peregrine, 2003).
The problem of wave-breaking in the ocean is multi-scale.Wave breaking must be represented in phase-average sea-state models which typically have the scale of (10-100 km) but ultimately these are trying to capture energy dissipation processes which occur at (1 μm-10 m) (Deike, 2022).Understanding the problem of wave breaking on a wave-by-wave basis (as opposed to a phased-averaged sea-state analysis) is frequently broken down into understanding when waves break and then understanding the subsequent evolution.Both parts of the problem are active areas of research.
Wave breaking criteria typically use the wave steepness, wave speed and water depth to provide a threshold above which waves would break.Various formulations for breaking criteria exist.Recent developments on this is the concept of breaking inception introduced in Derakhti et al. (2020), which describes the initiation of an irreversible process within the crest that leads to breaking and occurs prior to the stage of breaking onset.The breaking inception indicator proposed by Derakhti et al. (2020) is based on the diagnostic parameter  thresh used in the kinematic breaking criterion (Barthelemy et al., 2018).In Barthelemy et al. (2018), a threshold value  thresh ≈ 0.855 is suggested, and this experimental threshold value is verified in laboratory experiments and numerical studies for a variety of wave packet types.
The post-breaking behaviour is a complex multi-scale process.The breaking process generates turbulence, which is chaotic and complex.In simplified potential flow-based models this complexity can never (fully) be captured.However, we can attempt to model the key energy changes that occur when waves break.Liu et al. (2023) demonstrate that the cumulative dissipation effect of breaking itself is less chaotic and can be approximated.This suggests that prediction of the evolution of breaking waves (and thus wave forecasting) could be enhanced through improved physics-based models, data-driven methods, or other suitable techniques.
In addition to the strong non-linear process of wave breaking, weakly non-linear processes also play a role in water-wave evolution.One computationally fast model that captures this for narrowbanded https://doi.org/10.1016/j.oceaneng.2024.118009Received 7 February 2024; Received in revised form 22 April 2024; Accepted 22 April 2024 deep-water waves is the 'Modified non-linear Schödinger' (MNLS) equation (Dysthe, 1979) or variants of it (see also Dysthe and Trulsen, 2001).The MNLS equation models the evolution of the wave's complex envelope.Generally, these equations can predict the evolution of ocean waves well, although prediction of the largest waves remains challenging even for initially narrowbanded spectra due to the non-linear physics driving a local contraction of a wavegroup in space and which can also be thought of as a broadening of the spectrum (Adcock and Taylor, 2016).In random wave simulation there will always be some wave breaking, which is not captured in the basic MNLS equations.Thus, it is necessary to modify the standard equations in order to capture this process by adding a reduced-form (e.g., Trulsen and Dysthe (1990), Kato and Oikawa (1995), Hwung et al. (2011), Majda and Wang (2001) and Tian et al. (2012)) or data-driven breaking term (e.g., Liu et al. (2023)).Liu et al. (2023) showed that envelope-based reducedform breaking models can capture elements of wave breaking in the MNLS.but that this remains a challenging problem.
One approach to improving on envelope-based reduced-form breaking models is to use data-driven models and machine learning.In particular, we focus here on Artificial Neural Networks (ANN), highly efficient tools to capture the evolution of non-linear physics (Raissi et al., 2019;Kochkov et al., 2021).ANNs have been applied to study ocean waves in various ways (Qi and Majda, 2020;Eeltink et al., 2022).Attempts have been made to predict the evolution of significant wave height in the real ocean (Fan et al., 2020;Choi et al., 2020).Attempts have also been made to predict properties of the turbulence arising from wave breaking (Manucharyan et al., 2021) or the evolution of the interface in a two-phase system (Deo et al., 2023).Of particular importance here are applications to wave breaking.In Eeltink et al. (2022), the authors applied a blended model to predict wave breaking dissipation using the (temporal) MNLS equation (Dysthe, 1979;Shemer et al., 2010) combined with a long-short-term memory (LSTM) model (Hochreiter and Schmidhuber, 1997) to capture the effect of breaking dissipation.
In this study, we extend the work of Eeltink et al. (2022).We propose a novel blended model utilising a convolutional LSTM cell (ConvLSTM) (Shi et al., 2015) in conjunction with the spatial MNLS equation (Trulsen et al., 2000) to predict the evolution of focused wave-induced breaking dissipation in time.We aim at implementing a 'per-step' model to produce highly spatially resolved prediction in space.Synthetic focused wave data generated the envelope-based reduced-form breaking model of Liu et al. (2023) is applied to train this model.The long-term objective would be to train on highly spatially resolved Computational Fluid Dynamics (CFD) simulations or experimental measurements (if possible).We then analyse the performance of the trained model for focused waves, modulated plane waves and random waves.
The new approach has a number of advantages over (Eeltink et al., 2022) due to the introduction of the convolutional-recurrent architecture.We found that the ConvLSTM model learns breaking features locally, which enables the model to learn from a smaller set of representative samples before generalising these to arbitrary wave fields with flexible domain width and duration.The ConvLSTM model is able to predict the 'per-step' evolution of the wave breaking dissipation for highly spatially resolved wave field.Comparing this with the LSTM model, the behaviour of the ConvLSTM model is easier to interpret.In this study, we will show that our model is not limited to a single type of wave condition but can be effectively applied across various conditions even if not directly trained on these.
The present paper is part of a wider programme of work in which researchers are trying to use machine learning to capture multi-scale interactions with improved fidelity.In the present paper, we base our analysis on envelope-based models.The inability of these to capture wave breaking is one reason why their use has primarily been for academic research rather than as operational models.An reason for improving breaking characterisation is to make such a model more applicable to the real world and thus more useful operationally.
This paper is organised as follows.In Section 2 we introduce the method we use to generate the synthetic data.In Section 3 we introduce the model and training details.In Section 4 we evaluate the performance of the model using a test dataset not used for training.In Section 5 we analyse the performance of the for different types of waves, before summarising in the conclusions.

Methods
Our base model in the present paper is the Modified Nonlinear Schrödinger (MNLS) equation proposed by Trulsen et al. (2000); it solves the potential-flow equations and does not account for wave breaking.This model solves the linear part of the envelope evolution exactly and approximates the effects of weak nonlinearity under the assumption that the waves are narrowbanded.In the present work we assume the fluid is deep but not infinitely so, thus accounting for finite depth in the linear dispersion relation and in the return current but not in the coefficients of the non-linear terms which are only weakly influenced by depth for the depth range considered here (Barratt et al., 2021).

Envelope-based evolution equation
The velocity potential  and the surface elevation  are given by: where c.c. denotes the complex conjugate,  0 and  0 are the dimensional wavenumber and angular frequency of the carrier wave, (, , ) and (, ) are complex envelopes, and φ and η are the mean-flow potential and the wave-averaged free surface, respectively.The envelope (, ) is described by the MNLS equation (Trulsen et al., 2000) (given here in two dimensions or for unidirectional wave propagation).In Eqs.
(1) and ( 2), all variables are dimensional, and the primes as superscripts denoting their dimensionality are omitted here.It is more convenient to work with a non-dimensional system.The variables below are used as follows (Gramstad and Trulsen, 2011): where  is the characteristic steepness,  0 is the characteristic wave height,  0 is the length scale,  0 is the time scale,  0 is the nondimensional carrier wavenumber,  0 is the non-dimensional carrier angular frequency and  is the non-dimensional water depth.Variables with a prime indicating the corresponding dimensional variables.Variables without a prime are always non-dimensional (except for scaling variables  0 ,  0 and  0 ).Thus, we solve the following equations: where  * denotes the complex conjugate of , and  is the dissipation term that will be used to incorporate the reduced-form breaking models into the MNLS.In Fourier wavenumber space, the linear operator () can be written as (Trulsen et al., 2000): where  is the modulation wavenumber.The variables  and  follows the linear dispersion relationship: where  =  ′  2 0 ∕ 0 is the non-dimensional gravitational constant ( ′ = 9.81 m 2 s −1 ).
Eq. ( 4) is solved using a spectral method in space and timemarched using a fourth-order Runge-Kutta scheme.The complex envelope (,  = 0) is required as an initial condition for the solver.We either compute the envelope directly or reconstruct it from the Hilbert transform of the linearised surface elevation profile.In our simulation, we fix the steepness  = 0.35, so that  0 depends on the choice of carrier wavenumber.
The first-order surface elevation is obtained by shifting the complex envelope to the carrier frequency and taking the real part:

Reduced-form breaking model
We use the modified Kato & Oikawa (M-KO95 hereinafter) model developed in Liu et al. (2023).We also include in our comparison MNLS simulations with no breaking model (labelled 'NBM').The MNLS is able to integrate the equations even breaking would have occurred.
Kato and Oikawa (1995) (KO95) model breaking by modifying the return-current term.We use a spectral dissipation strength parameter β(), which is a function of wavenumber , to model the contribution to the dissipation term  on different wavenumbers.This model takes the form: where β() is obtained from fitting experimental measurements (see Liu et al., 2023) and the heaviside function  depends on whether a breaking criterion is satisfied.A global kinematic criterion with a fixed-duration active breaking interval is implemented.
The local kinematic breaking criterion, (∕ > 0.855 ± 0.05) from Barthelemy et al. (2018) and Derakhti et al. (2018), is based on the ratio of the horizontal fluid speed  and the crest speed .The former is obtained from the potential (cf. = ∕), which is reconstructed up to ( 3 ) according to (Carter et al., 2019): where B =   0  , we set  =  1 , and c.c. represents the complex conjugate.We ignore the effect of bound waves on the velocity as these are small.The breaking crest speed  is approximated by 0.8 0 ∕ 0 .The breaking criterion is local, which means it is applied only at those locations where ∕ > 0.855.
The global kinematic breaking criterion, (∕ > 0.855), will activate globally for a fixed time interval   whenever the kinematic breaking criterion is activated at any location, as detailed in Derakhti et al. (2018).We set the active breaking time interval   = 0.75 × (2∕(    )) (Derakhti et al., 2018), where we set   = 0.8 0 ∕ 0 and assume   =  local with the latter obtained from (Kurnia and Groesen, 2014): where  is a Hilbert transform.Choosing   when the kinematic breaking criterion (∕ > 0.855 ± 0.05) is satisfied at multiple locations, we take the maximum value and apply it globally.This treatment works for a scenario where only one breaking event is ever active at the same time.

Initial conditions for the MNLS simulations
We perform simulations for 3 types of wave conditions: (1) focused wave group (for training and testing purpose), (2) modulated plane wave (for testing purpose only), (3) random waves (for testing purpose only).We take the same sets of parameters from laboratory experiments on wave breaking for all types of waves (Eeltink et al., 2022) and give a brief review as follows.All variables in this section are dimensional with primes omitted.

Focused wave groups and random waves
We generate focused wave groups and random waves using the dimensional form of the JONSWAP spectrum (): where   is a scaling parameter to obtain appropriate initial wave amplitude  0 .All variables in are dimensional, but we omit primes for simplicity.For focused wave groups with JONSWAP spectrum, phases are assigned by computing a linear phase shift from a given focussing distance   based on linear theory.For random waves, the amplitudes are derived from the JONSWAP spectrum, and the phases are randomised by a fixed random number seed.As the final step, the surface elevation thus generated is made non-dimensional according to (3).The parameters for the focused wave group and random wave cases are drawn from the following ranges:  ∈ [2, 5],  0 ∈ [10, 60] mm,  0 ∈ [0.5, 1.25] Hz,   ∈ [10, 24] m for the focused wave only.

Modulated plane wave
We initialise our modulated plane wave case according to Eeltink et al. (2022) as where the upper and lower sideband frequencies are set according to  ± =  0 ± Δ with Δ the modulational frequency, and  is the phase shift of the sidebands compared to the carrier wave.Sideband wave numbers  ± are calculated using the linear dispersion relationship  2 ± =  ± tanh( ± ).The amplitudes of the carrier wave and the upper and lower sidebands are set according to   =  0 √   ,  + =  0 √  + , and , and  − = (1−  +)∕2, respectively, and   is the sideband fraction.All variables in Eq. ( 17) are dimensional, with all single quotes omitted.The parameters are drawn randomly from the following ranges:

Simulation settings
We perform non-dimensional MNLS simulations to generate synthetic wave breaking data.We use a 1D uniform mesh as our computational domain.The length of the computing domain spans [−40 0 , 40 0 ] in the direction of wave propagation, and the resolution in space is 0.15625 0 .The left-and right-hand sides of the domain are periodic boundaries.The duration of the simulation is 100 0 , and the resolution in time is  = 0.125 0 .We set the water depth to 0.8 m for all the cases in this study.

Neural network model
We introduce a novel approach to predict breaking dissipation locally in MNLS simulations, which is based on a convolution-recurrent neural network structure (Shi et al., 2015).In Fig. 1, we present the structure of the neural network model.It is composed of a ConvLSTM cell with tanh activation and a hard sigmoid for recurrent activation, a convolution layer, and an output layer, both with linear activation.We define the operators  ∶ Z → Z for the MNLS equation and the neural network model: where   MNLS indicates the operator for the MNLS equation ( 4) without a breaking model (NBM),   M−KO95 is the operator of the MNLS equation with the M-KO95 breaking model (hereafter MNLS+M-KO95) and   NN is the neural network operator. 0 is the initial condition and () is the solution to the initial value problem of Eq. ( 4) for  ∈ [0, ].() =   M−KO95 ( 0 ) −   MNLS ( 0 ) is the 'per-step difference' of the complex envelope and  is the time step.The 'per-step difference' is the change in the complex wave-envelope for each single time step. is the set of weights and biases of the neural network,  are the auxiliary states at time : where B is the approximation of  by Eq. ( 4), ,  are the hidden states and cell states of the neural network, respectively.Thus, we are seeking to model breaking as an infinitesimal-domain ML correction occurring over the real duration of the event.This contrasts with the previous work of Eeltink et al. (2022) where breaking was lumped.Periodic boundary conditions in space are embedded into the input and hidden-state tensor.
The model is trained with a complex envelope  as input and a perstep difference  as target.We use the Mean Squared Error (MSE) as loss function: where  batch ,   ,  are the batch size, grid number in space and grid number in time, respectively.We minimise the MSE loss according to Eq. ( 22).Since  is complex valued, Eq. ( 22) is identical in physics space and fourier space.Once the model is trained, we deploy the model for the prediction task.The algorithm for the simulation of the machine learning-corrected MNLS equation (hereafter MNLS+ML) is given by Algorithm 1.

4:
← B +   5: end for In Algorithm 1, we define several key variables that facilitate the execution of the MNLS+ML prediction: 1.   : This complex envelope evolves over iterations.

Focused wave cases
In this section, we explain the methodology used to generate the synthetic focused wave cases for the training phase of our model.Concurrently, we undertake a grid search to identify the optimal hyperparameters for our model.Subsequently, we analyse the performance of the model in predictive tasks.

Synthetic data generation
For training and evaluation purposes, we have synthesised a wave breaking dataset using the MNLS+M-KO95 model, which we designate as the 'ground truth'.This dataset comprises 1000 breaking focused wave cases, partitioned into 80% used for training, 15% used for validation, and used 5% for testing.We have chosen 1000 as a moderate number of cases and the model can generate acceptable predictions with fewer training cases, such as 100.We have also generated 50 modulated plane wave cases and random wave cases to investigate the generalisability of the model (see relevant discussion in Section 5).The parameter selection adheres to a uniform distribution across a range of parameters pertinent to ocean waves introduced in Section 2.3.In addition, we incorporate a classifier during the generation process to exclude cases characterised by an excessively high or low initial steepness.Cases with an inordinately high initial steepness often trigger wave breaking at the beginning of the simulation or could cause the solver to blow up due to the non-physical steepness.Conversely, cases with an excessively low initial steepness never trigger wave breaking throughout the simulation, reducing their utility for training.Both scenarios could potentially impede the learning rate.
Our training dataset comprises input and output fields.Initially, we perform the MNLS+M-KO95 simulation to obtain synthetic ground truth, denoted as  M−KO95 () =   M−KO95 ( 0 ), at discrete time steps  = 0, , 2, … , , where  represents the final time step's index and   =   • … •  ⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟  .Subsequently, we engage a 'per-step' NBM simulation at each time step, generating the per-step prediction  MNLS given by: Next, we compute the per-step difference between the ground truth and the estimated field, designated as () =  MNLS ()− M−KO95 (), for  = , 2, … , .For training the neural network model, the input data consists of a time series of complex envelope  MNLS (), at specific time steps  = , 2, … , .The output data, on the other hand, is defined by the change in the complex envelope, represented as the per-step difference (), at the same time steps.

Training history
Our model is trained for 50 epochs with an early stopping mechanism activated when there is no improvement in validation loss over 10 consecutive epochs.Moreover, we generate several models with different configurations to evaluate the impact of the choice of hyperparameters.The model is trained with data with time-steps of both 0.125 0 and 0.25 0 , and we do not find a clear link between performance and time resolution.The number and width of filters in the ConvLSTM layer turns out to be a critical hyperparameter that influences the model's size and performance.A grid search for these two variables indicates that our state-of-the-art (SOTA) model (which minimises the prediction loss for focused wave cases) employs 64 filters, each with a width of 9 (hereafter referred to as the w9f64 model) and comprises 160 898 trainable parameters.The training history of this SOTA model is illustrated in Table 1.Additionally, Table 1 shows the losses taken into account across all nine model configurations, revealing that the w9f64 model achieves the lowest prediction loss for the focused wave group cases.We also test a 'light' configurationthe w3f16 model -with a minimum of 16 filters of width 3 and a total of 5954 trainable parameters.This lighter configuration may be preferred when striking a balance between model accuracy and model size.The lighter configuration also seems to generalise to unseen wave types well, such as modulated plane waves and random waves.We also explored some other simple variants that enable stable prediction for long-term simulations.However, for the sake of simplicity, our study focusses primarily on the SOTA model, with a brief discussion of the other variants given in Appendix A.

Model analysis
In this section, we investigate the response of the neural network model to an arbitrarily selected case from the test dataset.Although we focus on this specific case, our conclusions generally apply to other scenarios.We specifically examine the cell states (Shi et al., 2015) of the ConvLSTM layer.In Fig. 2 we plot the spatial-temporal evolution of cell state for each channel (each filter in the ConvLSTM layer creates a channel, see Fig. 1).We find that some channels, especially the channels that maintain a high energy level throughout the duration, are activated over most of the simulation regardless of the state of the waves.Set against this, the other, 'weak', channels react to the input only during breaking.Interestingly, the cell states of some of the high-energy-level channels are blurred during breaking, which may be evidence that the model is switching between different states.This hints that the model has learnt to break the wave evolution problem into two processes -a classification process, which detects breaking, and a regression process, which occurs once breaking is detectedsimilar to how humans might tackle the breaking problem.The other evidence is obtained from the 'energy' of the cell state, defined as the mean squared sum over space:  , (, ) = (1∕  ) ∑    (,   , ) 2 of the ConvLSTM layer.In Fig. 3 we plot the evolution of the 'energy' of the cell state .The channel is sorted in descending order based on the cell energy at the last timestep.We find that the 'energies' of each channel are not equally distributed, and the 'energy' of some channels changes significantly during breaking.

Prediction for focused wave groups
We apply our neural network model to a test dataset of 50 focused wave cases unseen during the training phase.To illustrate its performance, we present a specific test case before assessing the error across the entire test dataset.This model has limited capabilities and inherits bias from the synthesised dataset, which makes it unsuitable to compare with laboratory experiments.However, we expect that a model trained with high-fidelity simulation data or experiment data can be compared with experiment measurements directly.
Fig. 4(a) compares the evolution of the complex envelope  in both the physical and wavenumber spaces for the NBM simulation, MNLS+ML simulation, and the ground truth.We note that the ground truth is synthetic data based on M-KO95 model, and for a thorough comparison among M-KO95 model, other reduced-order models, and laboratory data, can be found in Liu et al. (2023).Fig. 4(b) depicts the evolution of energy (‖‖ 2 ), peak frequency, and steepness for all three simulations, and the evolution of the MSE compared to the ground truth for NBM simulations and MNLS+ML simulations.Our observations in Fig. 4 suggest that the neural network model effectively predicts wave breaking dissipation and achieves a mean squared error (MSE) of the envelope approximately two orders of magnitude lower than that of the   In subplot (a), panels (i), (ii) and (iii) show the evolution of the absolute value of the complex envelope using different models, while panels (iv), (v), and (vi) show the evolution of the absolute value of the wavenumber spectrum of the complex envelope.Specifically, panel (i) and (iv) show the NBM result, panel (ii) and (v) show the MNLS+ML result, and panel (iii) and (vi) show the ground truth generated by the M-KO95 model.In subplot (b), panel (i) plots the mean squared error (MSE) of the absolute value of the complex envelope.Panels (ii), (iii), and (iv) compare the energy evolution  = ∫ ‖‖ 2 d normalised by its value at the first time step  0 , wavenumber spectral mean or centroid,  centre = ∫  B()d∕ ∫ B()d, and evolution of steepness.All panels share the same legend, which is placed on the top of subplot (b).In the legend, the green solid line represents the ground truth, the blue dashed line represents the MNLS+ML result or MSE of MNLS+ML compared to the ground truth, and the orange dotted line represents the NBM result or MSE of NBM result compared to the ground truth.
NBM simulation.The model also accurately reproduces the evolution of energy, peak frequency, and steepness.Fig. 5 presents the time series of the surface elevation profile and the frequency spectrum recorded at seven gauges within the domain.The recorded time series is shifted by the group speed calculated using linear dispersion, and subsequently cropped to a fixed duration  ∈ [0, 40].Our observations in Fig. 5 indicate that the predicted time series aligns closely with the ground truth in both the time and frequency domains, while the NBM solution quickly diverges from the ground truth.
Broadening our focus beyond a single instance, we assess the performance of the neural network model across all 50 focused wave group test cases.Fig. 6(a) presents the evolution of the mean squared error (MSE) for both the MNLS+ML simulations and the NBM simulations with respect to time.Further, Fig. 6(b) displays a histogram of the MSE at the final time step of the simulation.We observe that significant MSE errors arise from inaccurate predictions concerning both the onset and strength of the breaking event.

Modulated plane and random waves
To examine the model's generalisability beyond wave groups, we apply the trained model to predict breaking dissipation in modulated plane wave and random wave scenarios.We generate 50 modulated plane wave cases and 50 random wave cases for testing, the latter characterised by a JONSWAP spectrum.It is important to note that the performance of the M-KO95 model tends to deteriorate when applied to waves other than focused wave groups, owing to its non-local breaking dissipation.Consequently, given that the scenarios under consideration here are modulated plane waves and random waves, we anticipate observing larger discrepancies when compared to the true evolution as predicted by the M-KO95 model.
Figs. 7(a) and 7(b) plot the time series of surface elevation of a modulated plane wave and frequency spectrum.We find that the MNLS+ML simulation is much more accurate than the NBM simulations compared to ground truth.The error is, however, larger than for the focused wave cases presumably due to it being harder to localise the breaking event.Figs.7(c) and 7(d) plot the time series of surface elevation of  a random wave case and corresponding frequency spectrum and we are able to draw the same conclusion from the figures.We find that the neural network model is able to predict surface elevation with little discrepancy.Furthermore, Fig. 8 presents the evolution of the maximum breaking parameter, ∕, in the nonbreaking zone, averaged between test cases.We split nonbreaking and breaking zones in the MNLS+ML simulation by assessing whether the neural network output surpasses 3% of its total output range.Figs.8(a) and 8(b) compare the evolution of the mean breaking parameter, ∕, in MNLS+ML simulation with NBM and ground truth, for modulated plane wave and random waves respectively.The MNLS+ML and ground truth show similar ∕ evolution, while the NBM simulation exhibits unbounded growth.Error bars for the MNLS+ML series denote extremes of the breaking parameters in the test cases, with most maxima aligning or falling below the threshold value of 0.855 ± 0.05, validating that the MNLS+ML simulation adheres to the breaking criterion.As such, we can say that the ML architecture developed herein appears to work for other wave breaking scenarios and not just for isolated wave groups.Additional results are given in Appendix B.

Conclusion
Machine learning the complexities of wave breaking is a promising way of capturing the key physics of breaking in computationally fast, reduced-form models such as the MNLS equation.As such, an MNLS+ML model promises greater fidelity than alternative reducedform breaking models that have been implemented in the MNLS equation (see comparison of existing reduced-form breaking models in Liu et al. (2023)).This paper extends the framework proposed in Eeltink et al. (2022).The model proposed here achieves similar or superior performance compared to Eeltink et al. (2022), with a scalable structure that allows one to use from thousands to hundred thousands of parameters.The model utilises a convolutional-recurrent structure that enables it to support simulations with arbitrary length in space and time.
The model assumes that wave breaking dissipation is local in space.We train the model with focused wave cases only and show that the model is able to generalise to breaking in very different wave fields.The ability of a neural network model to learn from a small set of canonical samples and generalise to arbitrary conditions is promising.This reduces the potential difficulty in preparing a dataset for such models as both computational fluid dynamics (CFD) simulations and wave tank experiments for focused waves are cheaper than similar such simulations or experiments for modulated-plane-wave or random-wave datasets.
To improve the real-world accuracy of the model, integrating CFD data into the training process can be considered in future work, and the framework developed here is designed to enable this extension.This would enable the model to capture the complex dynamics of fluid flow and wave breaking.In addition, efforts should be made to establish a robust theoretical backbone for the model, incorporating insights from fluid dynamics, wave mechanics, and machine learning which the richer data from CFD might be able to provide.This would not only contribute to a deeper understanding of the underlying processes governing wave breaking and dissipation, but would also enable the development of more reliable and efficient models.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Figs. 9-11 provides a comprehensive review on the performance of the w3f16 model.We can see that for all types of waves, the model is able to predict the evolution of the surface elevation reasonably well.It can also forecast the evolution of energy, spectral centroid and steepness of the wave field.
Fig. 12 plots the structure of four simple variants of the base model we have explored in this study.In this figure, the models accept an input û and provide an residual output , where a shortcut is designed to build the final output by concatenating the input and the residual output in various ways.We find that after 200 epochs of training, the base model and variant V1 achieves a validation loss of 0.008, and variants V2, V3, V4 achieves a validation loss of 0.016.In addition, by utilising ReLU activation function, variants V2, V3 and V4 is unconditionally stable as these models only dissipates energy, which is suitable for long-term predictions.Notice that these unconditionally stable variants impose an assumption that the wave breaking only dissipates the energy and/or modifies the phases.

Appendix B. Additional results for plane wave and random waves
Fig. 13(a) demonstrates the evolution of the complex envelope  in both physical and wavenumber spaces across the NBM, MNLS+ML simulations, and the ground truth.Energy progression, peak frequency, and steepness for all three simulations, along with the MSE evolution for the NBM and MNLS+ML simulations relative to the ground truth, are depicted in Fig. 13(b).The neural network model, as inferred from Fig. 13, predicts wave breaking dissipation, outperforming the NBM simulation by reducing the envelope's MSE by two orders of magnitude.This model accurately reproduces the energy evolution, peak frequency, and steepness.A similar conclusion can be drawn from Figs. 14(a) and 14(b) for random waves, with the MSE, in this case, being 2.5 times lower than the NBM simulation.
Figs. 15(a) and 16(a) illustrates the temporal evolution of the mean squared error (MSE) for both the MNLS+ML simulation and the NBM simulation in the case of modulated plane waves and random waves, respectively.Additionally, Figs.15(b) and 16(b) presents a histogram of the MSE at the final time step of the simulation.Notably, the MSE shows a significant reduction of approximately an order of magnitude for modulated plane wave cases and 2.5 times for random cases.It is intriguing to observe that the model's performance in random cases correlates with the tail strength of the wave spectrum, which is governed by the bandwidth parameter .As the  decreases (i.e. the background spectrum becomes more broadbanded), the neural network prediction deteriorates rapidly.

Fig. 1 .
Fig. 1.Schematic of the structure and computing diagram of a neural network model composed of four layers in vertical: an input layer accepting a pair of real numbers (the magnitudes of the real and imaginary part of the complex envelope), a ConvLSTM layer, a convolution layer, and an output layer. is the timestep.The model accepts input B and model states  −1 ,  −1 from the previous timestep and predicts the change in   and model states   ,   between timesteps.
Y.Liu et al.

Fig. 3 .
Fig. 3.The evolution of the 'energy' (for each channel) of the cell state.The channels are sorted in descending order based on the energy level of the channel.

Fig. 4 .
Fig.4.The figures consists of two subplots.In subplot (a), panels (i), (ii) and (iii) show the evolution of the absolute value of the complex envelope using different models, while panels (iv), (v), and (vi) show the evolution of the absolute value of the wavenumber spectrum of the complex envelope.Specifically, panel (i) and (iv) show the NBM result, panel (ii) and (v) show the MNLS+ML result, and panel (iii) and (vi) show the ground truth generated by the M-KO95 model.In subplot (b), panel (i) plots the mean squared error (MSE) of the absolute value of the complex envelope.Panels (ii), (iii), and (iv) compare the energy evolution  = ∫ ‖‖ 2 d normalised by its value at the first time step  0 , wavenumber spectral mean or centroid,  centre = ∫  B()d∕ ∫ B()d, and evolution of steepness.All panels share the same legend, which is placed on the top of subplot (b).In the legend, the green solid line represents the ground truth, the blue dashed line represents the MNLS+ML result or MSE of MNLS+ML compared to the ground truth, and the orange dotted line represents the NBM result or MSE of NBM result compared to the ground truth.

Fig. 5 .
Fig. 5.The figure depicts profiles of a specific case from the test dataset at different gauge locations .The case is randomly selected from the test dataset.The solid green line is ground truth, the blue dashed line is the MNLS+ML simulation and the orange dotted line is the NBM simulation.Panel (a) plots the surface elevation  normalised by the characteristic wave amplitude  0 and panel (b) represents the frequency spectrum of normalised surface elevation, where  = ∕ 0 − 1 is the normalised modulational frequency.

Fig. 6 .
Fig. 6.This figure plots the MSE error distribution across the test dataset.Panel (a) shows the evolution of MSE of MNLS+ML and NBM simulations as a function of time, where the blue line with the shaded region corresponding to ±1 indicates the MSE of the MNLS+ML simulations and the orange line with the shaded region corresponding to ±1 indicates the MSE of the NBM simulations.Panel (b) illustrates the distribution of the MSE for the MNLS+ML and the NBM simulations at the final time step, where blue bars stand for MNLS+ML and orange bars for NBM.

Fig. 7 .
Fig. 7.The figure depicts (a) the surface elevation profile of a modulational plane wave case;(b) the corresponding frequency spectrum; (c) the surface elevation profile of a random wave case and (d) the corresponding frequency spectrum.The case was randomly selected from the test dataset.Solid green line is ground truth, blue dashed line is MNLS+ML and orange dotted line is NBM simulation.The time series of random wave case is not shifted by group speed.

Fig. 8 .
Fig. 8.The figure plots the evolution of maximum breaking parameter ∕ (Barthelemy et al., 2018) in nonbreaking zone averaged over test cases for (a) modulated plane wave and (b) random waves.Solid green line is ground truth, blue dashed line with errorbars is MNLS+ML and orange dotted line is NBM simulation.The errorbars indicate maxima and minima of maximum breaking parameter throughout the test cases.

Y
.Liu et al.

Fig. 11 .
Fig. 11.Various measurements taken from a randomly selected test cases of random wave, simulated by w3f16 model.The parameters for this case are  = 2.5328,  0 = 0.0383 m,  0 = 1.2343Hz and   = 1.2857.

Fig. 12 .
Fig. 12.This figure depicts four simple variants of the ConvLSTM model.

Fig. 13 .
Fig.13.The figure depicts various measurements taken from a randomly selected modulated plane wave case simulated by w9f64 model.In subplot (a), panels (i), (ii) and (iii) show the evolution of the absolute value of the complex envelope using different models, while panels (iv), (v), and (vi) show the evolution of the absolute value of the wavenumber spectrum.Specifically, panel (i) and (iv) show the NBM result, panel (ii) and (v) show the MNLS+ML result, and panel (iii) and (vi) show the ground truth generated by the M-KO95 model.In subplot (b), panel (i) plots the mean squared error (MSE) of the absolute value of the complex envelope.Panels (ii), (iii), and (iv) compare the energy evolution, wavenumber spectral mean or centroid, and evolution of steepness.All panels share the same legend, which is placed on the top of the subplot (b).In the legend, the green solid line represents the ground truth, the blue dashed line represents the MNLS+ML result or MSE of MNLS+ML compared to the ground truth, and the orange dotted line represents the NBM result or MSE of NBM result compared to the ground truth.

Fig. 14 .
Fig.14.The figure depicts various measurements taken from a randomly selected random wave case simulated by w9f64 model.In subplot (a), panels (i), (ii) and (iii) show the evolution of the absolute value of the complex envelope using different models, while panels (iv), (v), and (vi) show the evolution of the absolute value of the wavenumber spectrum.Specifically, panel (i) and (iv) show the NBM result, panel (ii) and (v) show the MNLS+ML result, and panel (iii) and (vi) show the ground truth generated by the M-KO95 model.In subplot (b), panel (i) plots the mean squared error (MSE) of the absolute value of the complex envelope.Panels (ii), (iii), and (iv) compare the energy evolution, wavenumber spectral mean or centroid, and evolution of steepness.All panels share the same legend, which is placed on the top of the subplot (b).In the legend, the green solid line represents the ground truth, the blue dashed line represents the MNLS+ML or MSE of MNLS+ML compared to the ground truth, and the orange dotted line represents the NBM result or MSE of NBM result compared to the ground truth.

Fig. 15 .
Fig. 15.This figure plots the MSE error distribution across the modulated plane wave test cases.Panel (a) shows the evolution of MSE of MNLS+ML prediction and NBM with respect to time where blue line with ±1 error bar indicates the MSE of MNLS+ML prediction and orange line with ±1 error bar indicates the MSE of NBM prediction.Panel (b) illustrates the MSE of MNLS+ML prediction and NBM at the final time step where blue bar stands for MNLS+ML and orange bar for NBM.

Fig. 16 .
Fig. 16.This figure plots the MSE error distribution across the random wave test cases.Panel (a) shows the evolution of MSE of MNLS+ML prediction and NBM with respect to time where blue line with ±1 error bar indicates the MSE of MNLS+ML prediction and orange line with ±1 error bar indicates the MSE of NBM prediction.Panel (b) illustrates the MSE of MNLS+ML prediction and NBM at the final time step where blue bar stands for MNLS+ML and orange bar for NBM.
: This variable serves as the ML correction.It is added to B to obtain the updated state   .This correction is part of the output derived from the   NN function.4.   and   : These auxiliary variables represent internal states that undergo updates at each iteration.They are associated with the   NN function, which updates these states simultaneously with   .
It is calculated by applying the MNLS equation, denoted by  MNLS , to the preceding state  −1 .Subsequently, a machine learning (ML) correction term   is incorporated into each iteration to adjust this state.2.B : This variable serves as aplaceholder, representing the tentative state of   before applying the ML correction.It is derived by applying the function   MNLS to the preceding state  −1 .Y. Liu et al.3.