Sunil Mohanlal

Jeffrey C Harris

Marissa L Yates

Stephan T Grilli

Unified depth-limited wave breaking detection and dissipation in fully nonlinear potential flow models

Keywords:

A new method is proposed for simulating the energy dissipation resulting from depth-limited wave breaking, in combination with a universal breaking onset criterion, in two-dimensional (2D) fully nonlinear potential flow (FNPF) models, based on a non-dimensional breaking strength parameter. Two different 2D-FNPF models are used, which solve the Laplace equation based on Chebyshev polynomial expansions or a boundary element method. In these models, impending breaking waves are detected in real time using a universal breaking onset criterion proposed in earlier work, based on the ratio of the horizontal particle velocity at the crest 𝑢, relative to the crest velocity 𝑐, 𝐵 = 𝑢∕𝑐 > 0.85. For these waves wave energy is dissipated locally using an absorbing surface pressure that is calibrated using an inverted hydraulic jump analog. This approach is first validated for periodic spilling breakers over plane beaches and bars, for which results are shown to be in good agreement with experimental data. Recasting this breaking dissipation model in terms of a non-dimensional breaking strength, the hydraulic jump analog is shown to provide results similar to those of a constant breaking strength model, and to yield good agreement for periodic plunging breakers as well. The same approach is then applied to irregular waves shoaling over a submerged bar, and shown to agree well with experimental data for the wave height, asymmetry, skewness, and kurtosis. Future work will extend this 2D breaker model to cases of three-dimensional (3D) breaking waves, simulated in existing 3D-FNPF models, in shallow or deep water conditions.

Introduction

Once generated by wind, ocean waves evolve, with complex kinematics and dynamics, as a result of nonlinear and dispersive effects, the effects of bathymetry, and dissipation from wave breaking and bottom friction, to name a few. Accurate simulations of this evolution are crucial for predicting phase-resolved surface wave properties in complex sea states, which govern wave interactions with fixed and floating objects, including offshore renewable energy systems, and surf zone parameters that drive nearshore currents and sediment processes, whose understanding and prediction are key to coastal management decisions.

Most wave processes that occur in complex sea states have already been simulated to some extent, based on equations representing the complete physics (i.e., derived from Navier-Stokes equations for single or multiple fluids). However, such simulations are highly computationally intensive and, hence, limited to small spatial and temporal scales. In contrast, operational models have been developed based on equations that simplify the wave physics, but nevertheless can simulate realistically many ocean wave processes over large areas and for long time periods. In such models, which are usually restricted to a specific wave regime (e.g., shallow or deep water conditions, small amplitude waves), important processes missing from the equations are parameterized in an ad-hoc manner, often on the basis of semi-empirical terms (e.g., breaking or bottom friction dissipation, the presence of structures). In this category are the standard phase-averaged ORCID(s):

wind wave models that are based on a spectral representation of the wave energy as a function of frequency and direction (e.g., STWAVE, [START_REF] Smith | STWAVE: Steady-state spectral wave model user's manual for STWAVE[END_REF]TOMAWAC, Benoit et al., 1997;WAVEWATCH III, Tolman, 2009).

However, in many engineering applications, phase-resolved wave properties are required in real time. Thus, there is a need for models in which wave elevation and kinematics are predicted over space and time, and in which the complex physics resulting from wave nonlinearity and fluid interactions with structures and the seafloor can be easily and accurately represented. Examples of problems requiring the use of such models include ship or ocean energy system seakeeping and motion control/optimization, coastal wave runup, and the prediction of extreme wave loads on ocean structures.

Phase-resolved models take many forms, depending on their domain of use. Ocean waves are often classified by the the depth 𝑑 to wavelength 𝐿 ratio. In deep water, for 𝑑 ≳ 𝐿∕2, where nonlinearity is usually weak but dispersive effects are important, one might apply models based on the linear mild-slope equation (MSE), such as REFDIF [START_REF] Kirby | A parabolic equation for the combined refraction-diffraction of Stokes waves by mildly varying topography[END_REF]. In shallow water, for 𝑑 ≲ 𝐿∕20, where nonlinearity and bottom effects dominate, and dispersive effects become less significant, models based on the Nonlinear Shallow Water (NSW) equations, which assume a uniform velocity over depth, could be an optimal choice (e.g., [START_REF] Stelling | An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation[END_REF][START_REF] Zijlema | Efficient computation of surf zone waves using the nonlinear shallow water equations with non-hydrostatic pressure[END_REF]. Lastly, in intermediate water depth, where both wave nonlinearity and dispersion are important, Boussinesq-type models [START_REF] Kirby | Boussinesq models and their application to coastal processes across a wide range of scales[END_REF] that feature both nonlinearity and dispersion to some extent, based on specifying a cutoff on higher-order terms representing these processes, are preferred provided that the horizontal velocity varies only moderately with depth such that it can be described by a polynomial approximation (e.g., [START_REF] Madsen | Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis[END_REF][START_REF] Agnon | A new approach to highorder Boussinesq models[END_REF][START_REF] Madsen | A new Boussinesq method for fully nonlinear waves from shallow to deep water[END_REF][START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking, and runup. I: 1[END_REF] and fully nonlinear Serre-Green-Naghdi-type models (e.g., [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: model development and analysis[END_REF][START_REF] Bonneton | A splitting approach for the fully nonlinear and weakly dispersive Green-Naghdi model[END_REF][START_REF] Shi | A highorder adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Zhao | On the steady solitary-wave solution of the Green-Naghdi equations of different levels[END_REF]. Boussinesq-type models are typically developed based on a perturbation expansion of the Fully Nonlinear Potential flow (FNPF) equations [START_REF] Kirby | Boussinesq models and their application to coastal processes across a wide range of scales[END_REF], with wave breaking, bottom friction, and horizontal vorticity effects represented by terms added to the equations to parameterize these physical processes [START_REF] Kennedy | Boussinesq modeling of wave transformation, breaking, and runup. I: 1[END_REF][START_REF] Shi | A highorder adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation[END_REF][START_REF] Kazolea | On wave breaking for Boussinesq-type models[END_REF]. These models have proved accurate in simulating laboratory experiments in which waves are generated in deep to intermediate water conditions and propagate into shallow water.

Models that directly solve the FNPF equations (e.g., [START_REF] Dold | Steep unsteady water waves: an efficient computational scheme[END_REF][START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF][START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF][START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF][START_REF] Grilli | A fully non-linear model for threedimensional overturning waves over an arbitrary bottom[END_REF][START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF][START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF][START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF][START_REF] Ducrozet | Applicability and limitations of highly non-linear potential flow solvers in the context of water waves[END_REF] are more computationally demanding than MSE, NSW, or Boussinesq-type models, but can accurately simulate waves in all water depth regimes up to wave breaking, since no assumptions are made about the wave nonlinearity or dispersion. FNPF models assume the flow to be irrotational and, hence, are governed by a Laplace's equation for the velocity potential, which makes them much more computationally efficient than full Navier-Stokes (NS) models. Among FNPF models, those based on the Higher-order Spectral (HOS) method [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF] are notably very efficient, but can typically only be applied to waves propagating in constant depth and up to a certain wave height, unless a modified form of the models is used, which makes them less efficient [START_REF] Ducrozet | Applicability and limitations of highly non-linear potential flow solvers in the context of water waves[END_REF]. FNPF models (except HOS-based models that require spatial periodicity) can also simulate wave interactions with structures and wave shoaling over an arbitrary bathymetry (e.g., slopes and/or bars), up to and into breaking/overturning (e.g., Grilli et al., 1994aGrilli et al., ,b, 1997[START_REF] Grilli | Closure of: Breaking criterion and characteristics for solitary waves on slopes[END_REF][START_REF] Grilli | Numerical modeling and experiments for solitary wave shoaling and breaking over a sloping beach[END_REF][START_REF] Grilli | Shoaling of periodic waves over barredbeaches in a fully nonlinear numerical wave tank[END_REF][START_REF] Guyenne | Numerical study of three-dimensional overturning waves in shallow water[END_REF][START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF][START_REF] Pomeau | Wave-breaking and generic singularities of nonlinear hyperbolic equations[END_REF]. However, FNPF models become unstable when waves begin to break and overturn, unless this process can be prevented by artificially specifying an energy dissipation for waves that are about to break (e.g., [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF] for a review).

As noted before, many NS models of wave breaking have been developed, using various numerical schemes and methods, including two-fluid models that represent the air and water (e.g., Guignard et al., 2001;[START_REF] Lachaume | Modeling of breaking and post-breaking waves on slopes by coupling of bem and vof methods[END_REF][START_REF] Abadie | Numerical simulation of waves generated by landslides using a multiple-fluid Navier-Stokes model[END_REF][START_REF] Banari | An efficient lattice Boltzmann multiphase model for 3D flows with large density ratios at high Reynolds numbers[END_REF][START_REF] Derakhti | Wave breaking in the surf zone and deep-water in a non-hydrostatic RANS model. part 1: Organized wave motions[END_REF], which can accurately simulate wave breaking either in direct NS simulations or based on standard turbulence models such as Large Eddy Simulation (LES; e.g., [START_REF] Harris | Large eddy simulation of sediment transport over rippled beds[END_REF]. However, the spatial resolution (and hence computational time) required to apply NS models currently restricts their use to small spatial and temporal scales and thus often to academic or idealized problems. Therefore, considering the large range of engineering applications involving strongly nonlinear and breaking waves, it is highly desirable to extend FNPF models to adequately model breaking waves and their related energy dissipation.

Although some earlier attempts to simulate the energy dissipation resulting from breaking waves in a FNPF model have been reported (e.g., [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wave-breaking onset[END_REF][START_REF] Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: wavebreaking energy dissipation[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF][START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], here a more general and accurate way of both detecting breaking onset in any conditions (i.e., wave types, bathymetry, and breaker types), including for nonlinear irregular wave trains, and simulating the corresponding energy dissipation in a more realistic manner are proposed. The present paper only describes twodimensional (2D) models featuring 2D breaking waves, but the extension to three-dimensions (3D) is in principle fairly straightforward and will be presented in future work.

In addition to being governed by a Laplace's equation, FNPF models use fully nonlinear kinematic and dynamic boundary conditions that are typically derived based on the formulations of [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] or [START_REF] Longuet-Higgins | The deformation of steep surface waves on water-I. A numerical method of computation[END_REF]. With the Eulerian framework of [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF], the free surface elevation is assumed single-valued, and thus waves can be modeled only until the instant that the free surface becomes vertical (e.g., [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF][START_REF] Craig | Numerical simulation of gravity waves[END_REF]Sulem, 1993, Bingham and[START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF]Yates andBenoit, 2015, Belibassakis and[START_REF] Belibassakis | A coupled-mode system with application to nonlinear water waves propagating in finite water depth and in variable bathymetry regions[END_REF]. In models that follow the Eulerian-Lagrangian framework of [START_REF] Longuet-Higgins | The deformation of steep surface waves on water-I. A numerical method of computation[END_REF], the free surface can be multivalued (e.g., for plunging breakers) and hence accurately simulated (in comparison to experiments) until the breaker jet impacts the free surface (e.g., in 2D, [START_REF] Dold | Steep unsteady water waves: an efficient computational scheme[END_REF][START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF][START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF]Grilli et al., 1997[START_REF] Grilli | Closure of: Breaking criterion and characteristics for solitary waves on slopes[END_REF][START_REF] Grilli | Numerical modeling and experiments for solitary wave shoaling and breaking over a sloping beach[END_REF]or, in 3D, Guyenne and Grilli, 2006;[START_REF] Fochesato | Numerical modeling of extreme rogue waves generated by directional energy focusing[END_REF].

In the present study, two 2D-FNPF models are modified and used to demonstrate the application of a novel combination of breaking onset/termination criteria and a breaker model: (i) the Eulerian model Misthyc [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], and the Eulerian-Lagrangian model, first proposed by [START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF], hereafter referred to as the "numerical wave tank" (NWT). Although the latter model can simulate all the cases presented here, some wave propagation cases can be more efficiently simulated by one FNPF model than the other. For instance, Misthyc, similar to other models of this type (e.g., [START_REF] Dommermuth | A high-order spectral method for the study of nonlinear gravity waves[END_REF][START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF][START_REF] Engsig-Karup | An efficient flexible-order model for 3D nonlinear water waves[END_REF][START_REF] Ducrozet | Applicability and limitations of highly non-linear potential flow solvers in the context of water waves[END_REF] that only apply to single-valued free surfaces, may be faster for certain domain sizes. In contrast, the higherorder Boundary Element Method (BEM) used in the NWT can simulate the exact geometry of a moving wavemaker, or complex bottom geometries, as well as multivalued free surfaces. Both types of models are equally able to simulate accurately highly nonlinear and dispersive waves propagating over arbitrary bathymetries. In both models, when waves are detected to approach breaking, an energy dissipation is explicitly specified to prevent wave breaking (e.g., instability or overturning) from occurring, which would interrupt the simulations. This is done by first identifying where and when in the computational domain impending breaking waves occur, using a breaking onset criterion, then applying a physically realistic energy dissipation onto these waves, commensurate with their parameters, and finally defining where and when this dissipation should cease to be applied, using a breaking termination criterion.

Wave breaking onset (or impending breaking) refers to the location in space and moment in time where and when some properties of an individual wave reach values indicating that the wave will soon begin to break and dissipate some of its energy through turbulence and viscosity. In a NS model of breaking waves, this process is included in the model equations and typically occurs automatically in the simulations. In models with reduced physics, such as FNPF, explicit breaking criteria must be defined to detect impending wave breaking. Such criteria usually depend on local wave properties such as their crest kinematics, steepness, surface slope, or curvature [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] reaching a specified threshold value. As the physics of wave breaking varies from deep to shallow water, until recently, different breaking criteria have been proposed to simulate different wave breaking regimes. In deep water, wave breaking is usually attributed to an excessive wave steepness, and referred to as steepness-limited breaking. In shallow water, wave breaking is usually induced by bathymetric effects, occurring during the shoaling process as waves propagate into shallow water, and is referred to as as depth-limited breaking. In the latter case, the type of wave breaking and the intensity of the related energy dissipation depend on the incident wave train and wave shoaling process, which is in general governed by the incident wave properties and bottom slope. Spilling breaking generally occurs over very mild slopes, plunging breaking over steeper slopes, and surging breaking over very steep slopes (see, e.g., Grilli et al., 1997 for solitary wave shoaling). Accordingly, many different definitions and criteria have been proposed for estimating breaking onset, that can be broadly classified [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] as geometric (e.g., [START_REF] Schäffer | A Boussinesq model for waves breaking in shallow water[END_REF], kinematic (e.g., [START_REF] Wei | A fully nonlinear Boussinesq model for surface waves. Part 1. Highly nonlinear unsteady waves[END_REF][START_REF] Kurnia | High order hamiltonian water wave models with wave-breaking mechanism[END_REF], or dynamic (e.g., [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] criteria. For depthlimited breaking, the geometric and kinematic criteria are often used (Grilli et al., 1997[START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF][START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], but they depend on an empirical constant that is a function of the bathymetry and incident wave conditions. Rather than relying on a collection of casespecific criteria, a universal breaking onset criterion is desirable, independent of empirical constants. In seminal work, [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] and [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] recently showed that there appears to be a universal breaking onset criterion for an evolving crest in the form of the ratio of the horizontal particle velocity at the crest to the wave (or crest) celerity, 𝐵 = 𝑢∕𝑐 reaching a threshold value 𝐵 𝑡ℎ = 0.85. This criterion indicates that, when the wave crest reaches this 𝐵 𝑡ℎ value, the wave will inevitably evolves towards breaking, although it does not necessarily start breaking at the limit of the threshold. This breaking onset criterion based on 𝐵 is used in the present work. It should be noted that [START_REF] Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wave-breaking onset[END_REF] and [START_REF] Seiffert | Simulation of breaking waves using the high-order spectral method with laboratory experiments: wavebreaking energy dissipation[END_REF] recently used this criterion in a HOS model, coupled with an energy dissipation model based on an eddy viscosity. Through laboratory experiments, they demonstrated the relevance and accuracy of this criterion for focused wave trains breaking over a flat bottom in intermediate water depths.

As energy dissipation resulting from wave breaking is complex and not yet fully understood, its simulation in models with reduced-physics, including those based on FNPF, has often relied on analogies with well-known dissipative phenomena, such as a hydraulic jump (HJ; e.g., [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF], and calibrating an empirical constant that may depend on the type of breaker. For example, a weak spilling breaker dissipates energy as white water rolls on the front face of the wave and, by moving in a frame of reference at the wave speed, this process resembles a bore or a HJ. Along this line, [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF] and [START_REF] Stive | Energy dissipation in waves breaking on gentle slopes[END_REF] compared the energy dissipated by a spilling breaker with that of a bore and estimated an empirical constant, 𝜇 ≃ 1.5, quantifying the ratio of energy dissipated by the breaking wave to that of an equivalent HJ. This analogy has been successfully used to simulate spilling breakers in FNPF models, in combination with various breaking criteria [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], Papoutsellis et al., 2019, Simon et al., 2019); it will also be applied in the present work. Although the HJ analog approach should be less accurate for plunging breakers, because breaking is more violent and the overturning surface no longer resembles a hydraulic jump, we shall see that it is still adequate in most cases. However, to simulate long-term irregular sea states, where both spilling and plunging wave breaking occurs, it is necessary to develop a breaker model that does not rely on test case dependent empirical constants.

Over the past few decades, various advances have been made in parameterizing wave breaking dissipation. [START_REF] Duncan | The breaking and non-breaking wave resistance of a two-dimensional hydrofoil[END_REF] conducted a series of experiments on steady breaking waves induced by fully submerged towed 2D hydrofoils, and quantified for these cases the non-dimensional breaking strength parameter 𝑏 (the wave breaking energy normalized by the fifth power of the wave celerity). Similarly, [START_REF] Phillips | Spectral and statistical properties of the equilibrium range in wind-generated gravity waves[END_REF], formulated spectral breaking strength as a function of wave speed using 𝑏 in deep water for irregular windgenerated waves in the ocean. Based on laboratory experiments, [START_REF] Romero | Spectral energy dissipation due to surface wave breaking[END_REF] followed up on this idea for deep water focused breaking waves, defining an empirical curve for 𝑏 as a function of the maximum surface slope of the focusing wave packet. Derakhti et al. (2018a) proposed an empirical relationship for parameterizing the breaking strength 𝑏 as a function of the time rate of change 𝑑𝐵∕𝑑𝑡 of the breaking onset parameter 𝐵 [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF], at breaking onset 𝐵 = 𝐵 𝑡ℎ , for deep or intermediate water 2D/3D focused waves. In this approach, by estimating the onset kinematics, the total energy dissipation of the resulting focused breaking can also be estimated. To utilize this in FNPF models, a universal breaking parameterization extending from shallow to deep water conditions is needed. In addition, the instantaneous energy dissipation needs to be explicitly specified. For these reasons, this work studies the parameterization of wave breaking in FNPF models, focusing on depth-limited breaking waves.

The paper is organized as follows. The formulations of the two FNPF models used are briefly described in Section 2. The breaking onset/termination criteria and the breaking dissipation models are described in Section 3, including the hydraulic jump model and a newer dynamic model proposed by Derakhti et al. (2018b). Applications of both models are presented in Section 4 for several regular and irregular depth-limited breaking wave cases, including some cases for which both types of breaking dissipation models are compared. Finally, the results are discussed and conclusions drawn in Sections 5 and 6.

FNPF models

FNPF models compute the irrotational motion of an incompressible and inviscid fluid, for which the fluid velocity u is represented by a scalar potential 𝜙, with u = ∇𝜙. For such flows, mass conservation becomes the Laplace equation for the potential,

∇ 2 𝜙 = 0 (1)
in the fluid domain Ω of boundary Γ.

For two-dimensional (2D) transient free surface flows in a vertical plane (𝑥, 𝑧), with a single-valued free surface elevation 𝜂(𝑥, 𝑡), the kinematic and dynamic free surface boundary conditions are,

𝜕𝜂 𝜕𝑡 = 𝜕𝜙 𝜕𝑧 - 𝜕𝜂 𝜕𝑥 𝜕𝜙 𝜕𝑥 (2) 𝜕𝜙 𝜕𝑡 = -𝑔𝜂 - 1 2 |∇𝜙| 2 - 𝑝 𝑎 𝜌 , (3)
on the boundary Γ 𝑓 and where 𝑔 denotes the gravitational acceleration, 𝜌 is the fluid density, and 𝑝 𝑎 is the free surface (atmospheric) pressure.

For typical 2D wave propagation problems, the fluid domain has an impermeable bottom boundary where a noflow, Neumann bottom boundary condition is specified as,

u ⋅ n = 𝜕𝜙 𝜕𝑛 = 0 (4)
on the boundary Γ 𝑏 , where n denotes the outward unit normal vector to the boundary. Typical conditions at the lateral boundaries of the 2D domain will be periodicity, a wave maker for generating waves, or an absorbing beach for dissipating waves. The details of these boundary conditions are provided in Section 4 for each specific application.

In the following subsections, the equations are briefly presented, and the numerical methods are summarized for the two existing FNPF models that use different versions of the free surface boundary conditions and different numerical methods for solving Laplace's Equation (Eq. 1).

Misthyc

The Misthyc FNPF model, developed by [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF], solves Laplace's Equation (Eq. 1) by mapping the potential 𝜙(𝑥, 𝑧, 𝑡) onto a boundary fitted vertical coordinate 𝑠 ∈ [-1, 1] and using a spectral approach to express 𝜙(𝑥, 𝑠, 𝑡) as a linear combination of Chebyshev polynomials (following [START_REF] Tian | A numerical model on the interaction between nearshore nonlinear waves and strong currents[END_REF]. At each time step, 𝜙(𝑥, 𝑠, 𝑡) is calculated by solving a system of 𝑁 𝑥 (𝑁 𝑇 + 1) linear equations, where 𝑁 𝑥 is the number of free surface nodes in the horizontal direction 𝑥 and 𝑁 𝑇 is the maximum order of the Chebyshev polynomials (here 𝑁 𝑇 = 7, following [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF].

Assuming single-valued free surface elevations 𝜂(𝑥, 𝑡), the free surface boundary condition Eqs. (2-3) are expressed following [START_REF] Zakharov | Stability of periodic waves of finite amplitude on the surface of a deep fluid[END_REF] as,

𝜕𝜂 𝜕𝑡 = w{ 1 + ( 𝜕𝜂 𝜕𝑥 ) 2 } - 𝜕𝜂 𝜕𝑥 𝜕 φ 𝜕𝑥 (5) 𝜕 φ 𝜕𝑡 = -𝑔𝜂 - 1 2 ( 𝜕 φ 𝜕𝑥 ) 2 + 1 2 w2 { 1 + ( 𝜕𝜂 𝜕𝑥 ) 2 } - 𝑝 𝑎 𝜌 , (6)
where φ(𝑥, 𝑡) = 𝜙(𝑥, 𝑧 = 𝜂, 𝑡) and w(𝑥, 𝑡) = 𝑤(𝑥, 𝑧 = 𝜂, 𝑡) = 𝜕 𝑠 𝜙(𝑥, 𝑠, 𝑡)| 𝑠=1 are the velocity potential and the vertical velocity on the free surface Γ 𝑓 , respectively (where subscripts indicate partial derivatives). Note that Eqs.

(5-6) can be transformed to express a relationship between 𝜂 and φ in the form of a so-called Dirichlet-Neumann (DtN) operator [START_REF] Craig | Numerical simulation of gravity waves[END_REF]. Following [START_REF] Bingham | On the accuracy of finite-difference solutions for nonlinear water waves[END_REF], [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] formulated the DtN problem corresponding to these equations. Assuming known 𝜂(𝑥, 𝑡) and φ(𝑥, 𝑡) values on Γ 𝑓 (𝑡), once w(𝑥, 𝑡) is computed from the solution of Laplace's equation, these quantities are advanced to time 𝑡 + Δ𝑡 by integrating Eqs. (5-6) with an explicit fourth-order Runge-Kutta scheme.

With this formulation in the Misthyc model, waves can only be simulated until the instant the free surface becomes vertical, at which time the model becomes unstable. To prevent this situation from occurring, impending wave breaking is detected using a breaking onset criterion, and a local damping is specified in the dynamic free surface boundary condition using an absorbing pressure 𝑝 𝑎 calibrated to simulate the wave breaking dissipation [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF][START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF]. Finally, waves are generated and absorbed at each end of the fluid domain by specifying relaxation zones that extend horizontally for ∼ 3𝐿, where 𝐿 is the dominant wavelength. Details are provided for specific cases in Section 4.

Numerical Wave Tank (NWT)

Longuet-Higgins and Cokelet (1976) first proposed a 2D-FNPF model to simulate overturning waves in a periodic domain with constant depth, until the instant the breaker jet impinges the free surface. The model solved Laplace's Eq. ( 1) at each time 𝑡, based on a complex potential Boundary Integral Equation (BIE) formulation (in a conformally mapped space) and integrated the kinematic and dynamic free surface boundary conditions Eqs. (2-3), expressed in a mixed Eulerian-Lagrangian frame of reference,

𝐷r 𝐷𝑡 = 𝜕r 𝜕𝑡 + (u ⋅ ∇)r = u = ∇𝜙 (7) 𝐷𝜙 𝐷𝑡 = -𝑔𝑧 + 1 2 |∇𝜙| 2 - 𝑝 𝑎 𝜌 , ( 8
)
by way of a predictor-corrector scheme, where r denotes the position vector on the free surface Γ 𝑓 . [START_REF] Dold | Steep unsteady water waves: an efficient computational scheme[END_REF] later proposed a more accurate time integration scheme for this model, based on an explicit Taylor series expansion of both r and 𝜙 on the free surface, which requires computing the successive material derivatives of Eqs. (7-8) and solving additional Laplace's equations for the corresponding time derivatives of the potential, up to a desired order. However, their model was still expressed in a conformally mapped domain and limited to space-periodic waves propagating in constant depth.

The 2D-FNPF model of [START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF]; [START_REF] Grilli | Corner problems and global accuracy in the boundary element solution of nonlinear wave flows[END_REF]; Grilli andSubramanya (1994, 1996); [START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF] used here is based on the same approach, but is formulated in the physical space, which allows for modeling wavemakers or other types of exact generation of fully nonlinear waves, such as from streamfunction wave theory, an absorbing beach at the far end of the domain, and an arbitrary bottom bathymetry. Laplace's Eq. ( 1) is solved based on a BIE derived from Green's second identity,

𝛼(𝐱 𝑖 )𝜙(𝐱 𝑖 ) = ∫ Γ { 𝜕𝜙 𝜕𝑛 (𝐱)𝐺(𝐱 -𝐱 𝑖 ) -𝜙(𝐱) 𝜕𝐺 𝜕𝑛 (𝐱 -𝐱 𝑖 ) } dΓ, ( 9
)
where 𝜙 is the velocity potential on the boundary Γ, 𝛼 is the interior angle made by the boundary at point 𝐱 𝑖 and 𝐺(𝐱, 𝐱 𝑖 ) = -(1∕2𝜋) ln 𝑟 𝑖 is the 2D free space Green's function (with 𝑟 𝑖 = |𝐱 -𝐱 𝑖 |). In the model, this equation is discretized by various types of higher-order boundary elements and both regular and singular, as well as quasisingular (occurring when two parts of the boundary are close to each other, e.g., in the tip of breaker jets) integrals are computed by very accurate methods (Grilli andSubramanya, 1994, 1996). Additionally, extended compatibility conditions of the solution on both sides of the boundary are specified at corners in the domain (e.g., between a wave maker and the free surface, [START_REF] Grilli | Corner problems and global accuracy in the boundary element solution of nonlinear wave flows[END_REF][START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF].

Assuming known 𝜂(𝑥, 𝑡) and 𝜙(𝑥, 𝑡) values on Γ 𝑓 (𝑡), once 𝜕 𝑛 𝜙 is computed from the solution of Eq. ( 9), these quantities are advanced to time 𝑡+Δ𝑡, by integrating Eqs. (7-8), as in [START_REF] Dold | Steep unsteady water waves: an efficient computational scheme[END_REF], based on explicit Taylor series expansions of both r and 𝜙, limited here to secondorder, hence requiring to solve an additional Laplace's equation for 𝜕 𝑡 𝜙. This is done using a similar BIE to Eq. ( 9) in the same discretized computational domain, for a modest additional computational effort.

To generate waves, the NWT can simulate the motion of a flap or piston wavemaker, or the generation of numerically exact streamfunction waves on a lateral boundary [START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF]. For wave absorption, an absorbing beach can be specified at the far end of the domain, combining an absorbing free surface pressure for high-frequency waves and an absorbing lateral piston wavemaker for lowfrequency waves [START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF].

Wave breaking model

Wave breaking in a FNPF model requires three steps: (i) a breaking onset criterion, which allows identifying where on the free surface and when a wave has reached a threshold beyond which breaking is inevitable, and indicating where on the free surface and when the application of an energy absorption in the model should start; (ii) based on wave parameters, a method to quantify the energy dissipation rate that should be specified in the model to damp waves that are identified to be breaking in (i); and (iii) a breaking termination criterion, which indicates where on the free surface and when the application of the energy dissipation should end.

Breaking onset criterion

As discussed in the introduction, in this work, breaking onset is detected based on the universal criterion proposed by [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] for deep and intermediate water depth breaking, and validated by [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] for shallow water breaking of any type (i.e., spilling, plunging, or surging). These studies showed, for a wide variety of conditions leading to breaking such as energy focusing or effects of bathymetry, that a steepening wave whose ratio of the horizontal particle velocity at the crest 𝑢 to the wave (or crest) celerity 𝑐, noted 𝐵 = 𝑢∕𝑐, exceeds the threshold value 𝐵 𝑡ℎ = 0.85, will eventually break; and waves for which 𝐵 < 𝐵 𝑡ℎ will not break. It should be emphasized that wave breaking does not start at this threshold, but later on when 𝐵 ≃ 1. Instead, the 𝐵 criterion predicts where and when a wave crest passes a point of no return, beyond which it will eventually break. [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] showed that, in shallow water, this occurs about 0.2𝑇 in average after breaking onset, where 𝑇 is the mean wave period in the sea state. In the FNPF models, it takes time for the absorbing pressure that is used to damp impending breaking waves (see next subsection) to extract enough energy from the wave crest to prevent breaking. Therefore, the energy dissipation in the models is applied when a wave reaches the threshold 𝐵 = 𝐵 𝑡ℎ .

One challenging aspect of applying this kinematic breaking criterion is that it is local and requires first that all wave crests be identified and tracked at all times over the entire free surface Γ 𝑓 , and second that the crest celerity, 𝑐 = 𝑑𝑥 𝑐 ∕𝑑𝑡 is computed accurately (where 𝑥 𝑐 (𝑡) denotes a given wave crest location as a function of time). Accordingly, in the models, similar to [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF], [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], and [START_REF] Stansell | Experimental investigation of wave breaking criteria based on wave phase speeds[END_REF], the instantaneous celerity 𝑐(𝑡) of individual tracked wave crests is computed by fitting a moving polynomial to 𝑥 𝑐 (𝑡) over a

Breaking dissipation

Absorbing surface pressure

For waves that have been identified as evolving to breaking by the breaking onset criterion, the energy dissipation is specified, as in earlier work [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], using an absorbing (or damping) surface pressure 𝑝 𝑎 in the dynamic free surface boundary condition, Eq. ( 6) or (8) for the Misthyc and NWT models, respectively. This pressure is applied spatially across part of the back and front of the breaking wave crest, and it is defined as being proportional to the normal velocity of water particles at the free surface (Fig. 1),

𝑝 𝑎 (𝑥, 𝑡) = 𝜈 𝑎 (𝑡) 𝑆(𝑥) 𝜕𝜙 𝜕𝑛 (𝑥, 𝑡) ∀𝑥 ∈ (𝑥 𝑙 , 𝑥 𝑟 ) (10)
where 𝑆(𝑥) is a non-dimensional shape function (defined below) and 𝜈 𝑎 (𝑡) is an absorption function, with the dimension of 𝜌 𝑐 (i.e., units of kg m -2 s -1 ), defined such that the rate of work produced by the pressure against the wave motion, 𝑝 𝑎 𝜕 𝑛 𝜙, integrated over the selected segment of the breaking wave surface, matches an expected rate of energy dissipation per unit length of crest, Π 𝑏 for the wave, as will be described in the following subsection.

For single valued surface elevations Eq. ( 10) yields,

𝜈 𝑎 (𝑡) = Π 𝑏 ∫ 𝑥 𝑟 𝑥 𝑙 𝑆(𝑥)(𝜕 𝑛 𝜙) 2 √ 1 + (𝜕 𝑥 𝜂) 2 𝑑𝑥 . ( 11
)
in which the square root is the Jacobian, 𝑑Γ∕𝑑𝑥.

The limits of integration in Eq. ( 11) are defined across each breaking wave crest, with (𝑥 𝑙 , 𝑥 𝑟 ) located near the troughs located behind and ahead of the breaking wave crest, respectively, such that |𝜕 𝑛 𝜙| < 𝜀 |𝜕 𝑛 𝜙| 𝑚𝑎𝑥 , where 𝜀 ≪ 1 (here 𝜀 = 10 -4 ). As shown in Fig. 1, 𝑆(𝑥) is a smooth function that varies from 0 to 1 over the breaking region, with a ramp to ensure smooth transitions of 𝑝 𝑎 between zero in non-breaking regions and its calculated value in breaking regions [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF],

𝑆(𝑥) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0, 𝑥 ≤ 𝑥 𝑙 cos ( 𝜋 2 𝑥-𝑥 𝑙1 𝑥 𝑙 -𝑥 𝑙1 ) , 𝑥 𝑙 ≤ 𝑥 ≤ 𝑥 𝑙1 1, 𝑥 𝑙1 ≤ 𝑥 ≤ 𝑥 𝑟1 cos ( 𝜋 2 𝑥-𝑥 𝑟1 𝑥 𝑟 -𝑥 𝑟1 ) , 𝑥 𝑟1 ≤ 𝑥 ≤ 𝑥 𝑟 0, 𝑥 ≥ 𝑥 𝑟 ( 12
)
where 𝑥 𝑙1 = 𝑥 𝑙 + 𝛼(𝑥 𝑟 -𝑥 𝑙 ) and 𝑥 𝑟1 = 𝑥 𝑟 -𝛼(𝑥 𝑟 -𝑥 𝑙 ), 𝛼 = 0.1.

Rate of energy dissipation

Two different parameterizations of the rate of energy dissipation in breaking waves, Π 𝑏 , will be considered in this work. The first one, based on the hydraulic jump analogy introduced in earlier work [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF][START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF][START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], will be the default method used in both models, and the second one, based on the time rate of change of the breaking onset criterion, recently proposed by Derakhti et al. (2018b), will be used for comparison.

Hydraulic jump analogy:

In earlier work, [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF], [START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF], and [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF] parameterized Π 𝑏 based on the classical energy dissipation of an hydraulic jump (HJ), i.e.,

Π ℎ = 𝜌𝑔𝑐 𝑑 𝐻 3 4ℎ 𝑐 ℎ 𝑡 with Π 𝑏 = 𝜇Π ℎ , ( 13
)
where 𝑐 is the wave phase speed (or crest celerity), 𝑑 the undisturbed water depth below the point of maximum front slope, 𝐻 the wave height (measured trough to crest), ℎ 𝑐 the total depth below wave crest, and ℎ 𝑡 the total depth below wave trough (Fig. 1). As shown in [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF] and [START_REF] Stive | Energy dissipation in waves breaking on gentle slopes[END_REF], the equation for Π ℎ can be obtained easily from a control volume approach by deriving equations for mass, momentum, and energy conservation assuming periodic waves with a uniform velocity over depth, hydrostatic pressure, and negligible bottom friction. Based on laboratory experiments for spilling breakers propagating over mild slopes, [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF] proposed that 𝜇 = 1.5. Note, for symmetric linear waves breaking over mild slopes, 𝐻 = 2𝑎, with 𝑎 the wave amplitude, 𝑐 ≃ 𝑐 𝓁 = √ 𝑔𝑑, ℎ 𝑡 ≃ 𝑑 -𝑎 and ℎ 𝑐 ≃ 𝑑 + 𝑎, and Eq. ( 13) transforms into,

Π 𝑏 = 𝜇 ( 𝜖 𝑏 ) 2𝐹 3 1 -𝐹 2 (14) with 𝜖 = 𝑏𝜌𝑔 -1 𝑐 5 , ( 15
)
where 𝐹 = 𝑔𝑎∕𝑐 2 𝓁 is the wave Froude number [START_REF] Kirby | Discussion of 'Note on a nonlinearity parameter of surface waves' by S. Beji[END_REF], which for long breaking waves in shallow water reduces to 𝐹 ≃ 𝑎∕𝑑, 𝜖 is the breaking energy dissipation rate per unit of wave crest width proposed by [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF][START_REF] Duncan | The breaking and non-breaking wave resistance of a two-dimensional hydrofoil[END_REF], based on measurements of steady spilling breakers in deep-water, and 𝑏 is a non-dimensional breaking strength parameter that was parameterized as a function of the hydrofoil characteristics and submergence depth used in these experiments. To obtain Π 𝑏 = 𝜖 requires that 𝑏 ≃ 2𝜇𝐹 3 in Eq. ( 14), assuming 𝐹 2 ≪ 1. Based on Svendsen's parameterization, for small 𝐹 , the HJ breaking strength parameter can thus be expressed as

𝑏 ≃ 2𝜇𝐹 3 = 3𝐹 3 .
Making various scaling arguments of the turbulent energy dissipation for unsteady breakers in deep or intermediate water, [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF] showed that 𝑏 ∝ 𝐹 2.5 , while more recently, using different scaling arguments, [START_REF] Mostert | Inertial energy dissipation in shallow-water breaking waves[END_REF] proposed that 𝑏 ∝ 𝐹 3.5 for unsteady shallow water breakers. Thus, the above formulation of 𝑏 falls in between these independent results, without any clear guidance on which result is most realistic. Note that, consistent with the analysis of [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF], [START_REF] Romero | Spectral energy dissipation due to surface wave breaking[END_REF] and Derakhti et al. (2018b) assumed that 𝑏 ∝ 𝐹 2.5 . It is outside the scope of this paper to attempt to reconcile these conflicting results, which will be the object of other studies [START_REF] Derakhti | Predicting the rate of energy dissipation in breaking gravity water waves from deep to shallow water[END_REF].

The default parameterization of energy dissipation in both FNPF models used in the present applications will thus be based on Eq. ( 13), with 𝜇 = 1.5. It should be noted that, similar to the breaking onset criterion used to detect impending breaking in the models, computing Π ℎ requires identifying individual wave crests and troughs and accurately computing the crest celerity. The same method used to compute the onset criterion is thus used to compute the parameters required for estimating the wave energy dissipation.

Parameterization based on time rate of breaking onset criterion: By performing numerical simulations of focused waves, [START_REF] Romero | Spectral energy dissipation due to surface wave breaking[END_REF] Derakhti et al. (cf., Fig. A2 in Derakhti et al., 2020). In some of the applications in Section 4, 𝑏 will be computed with Eq. ( 16) and compared to results based on the HJ analog. To do so, the time rate of change of 𝐵 at the threshold will be calculated by applying a linear fit to the 𝐵 values computed in the interval [0.82-0.85]. Note that 𝑏, which quantifies the average energy dissipation rate during a breaking event based on 𝛾, computed at breaking onset, is not the time average of 𝑏(𝑡), which quantifies the instantaneous energy dissipation rate based on instantaneous wave parameters.

As discussed before, [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] showed that the 𝐵 criterion also applies to shallow water breaking waves and, hence, is universal. Calculating the energy dissipation rate in their model for shallow water breaking waves, Derakhti et al. (2018b) confirmed the parameterization of 𝑏 in Eq. ( 16) for 𝛾 < 1.11, with the breaking strength parameter reaching an upper bound 𝑏 𝑚𝑎𝑥 = 0.02 for larger 𝛾 values, although this parameterization, particularly for large 𝛾, was revised by [START_REF] Derakhti | Predicting the rate of energy dissipation in breaking gravity water waves from deep to shallow water[END_REF].

Given Π 𝑏 = 𝜖, the second parameterization of energy dissipation considered in both FNPF models in the present applications will be based on the expression of 𝜖 in Eq. ( 15), with 𝑏 obtained from Eq. ( 16) for 𝛾 ≤ 1.3, and 𝑏 = 𝑏 𝑚𝑎𝑥 for 𝛾 ≥ 1.3. As before, most wave crest parameters required to compute Π 𝑏 are similar to those used to compute 𝐵 and are already available in the models.

Breaking termination criterion

Breaking termination is also an important factor to extract accurately the appropriate amount of energy from breaking waves in the models. However, unlike in actual waves, as would for instance be simulated in a NS-VOF model [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF], the value of 𝐵 does not grow much beyond or remain above 𝐵 𝑡ℎ in the FNPF model once dissipation is applied and a different method is required to detect breaking termination. To be consistent with the onset criterion, the termination criterion is based on 𝐵 reaching a value 𝐵 off that is lower than the onset threshold, and needs to be calibrated based on benchmark data. For the applications considered here, the optimal 𝐵 off value appeared to be problem dependent. Specific values and their implications are discussed in the section 4.2.

Applications

Comparison of the two formulations of breaking dissipation

Before presenting detailed applications using the same breaking onset/termination criteria, it is of interest to analyze and compare values of 𝐵 and 𝑑𝐵∕𝑑𝑡 computed near and at the breaking onset for various cases, as well as the breaking strength parameter 𝑏 and/or dissipation rate Π 𝑏 resulting from the two parameterizations of the rate of energy dissipation discussed before. These are the experimentally validated depth-limited HJ spilling breaker model [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], whose energy dissipation rate is given by Eq. ( 13) and for which 𝑏 can readily be obtained based on the expression of 𝜖 in Eq. ( 15), assuming Π 𝑏 = 𝜖, and the newer dissipation rate based on 𝑑𝐵∕𝑑𝑡 and 𝛾 (Derakhti et al., 2018b), with 𝑏 given by Eq. ( 16).

In the following, breaking wave parameters calculated based on the HJ breaking model, with Misthyc and the NWT, are compared in simulations of laboratory experiments for periodic: (i) shoaling and plunging breaking waves propagating over a bar, from [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF] In each test case, the HJ breaking dissipation model is used, and the parameter 𝑏 is calculated using the expression of 𝜖 in Eq. ( 15), based on the average power dissipated in the model over the breaking onset/termination range, 𝜖 = ⟨Π 𝑏 ⟩, through the application of the absorbing pressure 𝑝 𝑎 based on Eqs. (10-12), with the instantaneous Π 𝑏 (𝑡) given by Eq. ( 13). For the second parameterization of dissipation, the corresponding 𝑏 values are found using Eq. ( 16), based on the wave parameters computed at breaking onset. Recall also that the breaking strength shown here is based on an average energy dissipation rate, ⟨𝜖⟩, which will have a different value depending on the duration of active breaking, which in both models could be adjusted without affecting significantly the resultant wave characteristics. Finally, in this paper potential flow is assumed, while wave breaking transfers energy to non-potential components of the flow (e.g., vorticity) that may not be dissipated and would appear as an over-prediction of the energy dissipation.

Based on the 𝑏 values simulated in these applications, an instantaneous dissipation rate Π 𝑏 = 𝜖 defined with the expression of 𝜖 in Eq. ( 15), based on a constant breaking strength 𝑏 = 0.05, was also tested. Since the wave celerity used in the definition of 𝜖 is a function of space and time, 𝑐(𝑥, 𝑡), this parameterization also provides a time-varying dissipation rate. For instance, for the BB-regular case, Fig. 4 shows that this constant breaking strength yields an instantaneous dissipation rate that agrees well with that calculated using Eq. ( 13) for the HJ parameterization. Therefore, for the BB-regular case, the constant strength approach would likely be accurate enough, and Fig. 3 shows that for the wider range of cases that will be considered hereafter, the average 𝑏 value computed in the model in each case is on the same order as 𝑏 = 0.05.

Experimental validations

In this subsection, the results of numerical simulations with one or both FNPF models, using the energy dissipation rate of the HJ model, the constant breaking strength 𝑏 = 0.05, or both, are compared with laboratory experiments performed for five standard benchmark cases. The five test cases are the periodic spilling breaker experiments over a plane slope of [START_REF] Hansen | Regular Waves in Shoaling Water Experimental Data[END_REF], and [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF], the regular and irregular plunging breakers 13), or (red) the Mis-b005 Eq. ( 15) parameterization with 𝑏 = 0.05. over a bar of [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF], and the irregular breaking waves over a bar of [START_REF] Adytia | Numerical simulation of breaking regular and irregular wave propagation above a sloping bottom[END_REF]. In the models, fully nonlinear periodic waves are generated based on streamfunction wave theory in the generation zone for Misthyc [START_REF] Benoit | Near-bottom kinematics of shoaling and breaking waves: Experimental investigation and numerical prediction[END_REF], and using an exact wavemaking (particle curtain) boundary [START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF] for the NWT. Irregular waves are simulated in the model using a flap wave maker boundary (e.g., [START_REF] Grilli | Numerical generation and absorption of fully nonlinear periodic waves[END_REF].

The model discretizations in space and time are specified such that the Courant number, CFL ≈ 1.0 in Misthyc (as suggested by [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and CFL ≈ 0.45 in the NWT (found to be optimal by [START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF]. Using the optimal CFL number for each model, a refined spatial discretization is specified in each case to ensure high numerical accuracy (see [START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF] for guidance). Fig. 5 shows examples of the instantaneous free surface elevation computed for the HS test case for a range of spatial discretizations after simulations reach a quasisteady state in the models. These are defined based on the initial spatial discretization on the free surface Δ𝑥, where 𝐿∕Δ𝑥 = 35, 50, or 70, and 𝐿 is the incident wave length. The simulation results show that free surface elevations are nearly identical in deeper water for both models over the range of tested discretizations. Small differences can only be seen in shallower depths, particularly near the wave crests. Based on these results, in all the simulations discussed hereafter, the spatial discretization was prescribed such that 𝐿∕Δ𝑥 > 50 (note, for irregular waves, 𝐿 denotes the dominant wavelength). Table 2 summarizes the spatio-temporal parameters used in the simulations with the Misthyc and NWT models for the five benchmark cases detailed in the following sections.

Periodic spilling breakers on a slope -Hansen and Svendsen (1979)

Hansen and Svendsen (1979) (referred to as HS) performed experiments for periodic waves shoaling and spilling breaking waves propagating over a mild slope. The wave tank had a constant initial depth of ℎ 0 = 0.36 m, from the wavemaker up to 𝑥 = 14.78 m, the toe of the 1∕34.26 slope. Regular waves with an initial height 𝐻 0 = 0.095 m, period 𝑇 = 1 s, and incident wavelength 𝐿 0 = 1.43 m were generated at the wavemaker (note, these waves were generated in intermediate waver conditions, with ℎ 0 ∕𝐿 0 = 0.252). This benchmark was simulated with both Misthyc and the NWT, using the energy dissipation rates from the HJ and Mis-b005 or NWT-b005 models, respectively (see Table 2 for the numerical parameters used). Since wave runup was not considered in either model, an absorbing beach was modeled for 𝑥 ≥ 25.5 m with a deepening bathymetry (for 𝑥 ∈ [25.5 -27] to induce deshoaling, which aids the absorption of waves) followed by constant depth in the absorption zone (see Fig. 6a and Fig. 1 in Grilli et al., 2020 for details). In these experiments, the breakers reach the shoreline, so 𝐵 off = 0 is used as the breaking termination criterion. Preliminary tests with larger values appeared to cause wave reformation that was not observed in the experiments. In the numerical models, breaking onset with 𝐵 = 𝐵 𝑡ℎ = 0.85 occurs at 𝑥 𝑏 ≃ 22.2 m, as compared to 𝑥 𝑏 ≃ 22.5 m in the experiments. The model results were averaged over 5 successive wave periods after the simulations reached a quasi-steady state Fig. 6b shows (in both the experiments and all numerical simulations) the wave celerity normalized by the deep water linear celerity 𝑐 0 = 𝑔𝑇 ∕(2𝜋) gradually decreases in the shoaling region. In general, the simulation results agree well with the experimental measurements, but less so in breaking region (𝑥 > 𝑥 𝑏 ; where experimental data is very noisy) due to unsteady variations in calculations of the derivative of the wave crest displacement. For 𝑥 < 𝑥 𝑏 , 𝑐∕𝑐 0 is slightly larger than in experiments, which is consistent with the results of [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF]. Fig. 6c shows the wave height normalized by 𝐻 0 gradually increases over the shoaling region then rapidly decreases beyond breaking onset and, in all cases, the simulation results agree well with experimental measurements.

Considering the two parameterizations of energy dissipation in Misthyc, the simulation results show only small differences throughout the simulations, and these differences, as well as those with the NWT results, mostly occur at the far end of the tank for 𝑥 > 23.5. The differences likely result from the different numerical methods and discretizations, as well as the regridding used in the NWT, where the Eulerian-Lagrangian approach causes grid points to cluster around the breaking crests. To limit this clustering, regridding is periodically calculated in the model using cubic shape functions [START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF], which may slightly affect the application of the breaker model around the breaking crests.

Note, when using the NWT with the HJ model and a geometric breaking criterion (front slope 𝛽 𝑚𝑎𝑥 = 37 𝑜 ) for this case, [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF] predicted breaking onset slightly

Periodic spilling breakers on a slope -Ting and Kirby (1994)

Experiments of shoaling and spilling breaking waves propagating over a plane 1/35 slope, similar to those by HS in the previous section, were performed by [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF] (referred to as TK). Periodic waves of height 𝐻 0 = 0.125 m, period 𝑇 = 2.0 s, and wavelength 𝐿 0 = 3.85 m, were generated by a piston wavemaker and propagated in a tank of constant depth ℎ 0 = 0.4 m to 𝑥 = 3.85 m, the toe of the slope. Similar to the previous test case, in the models, an absorbing beach was specified for 𝑥 > 15.5 m, starting in a water depth ℎ = 0.067 m and gradually deepening to ℎ = 0.29 m at 𝑥 = 18 m. As in the previous case, 𝐵 off = 0. Fig. 7 shows the experimental setup along with the the locations of the 12 wave gauges that measured the free surface position.

For this test case, the models detected the onset of wave breaking at 𝑥 𝑏 ≃ 10. using the 𝐵 = 0.85 criterion, both models predict more accurately the wave elevation at this gauge. Beyond breaking (for 𝑥 > 𝑥 𝑏 ), results from Misthyc using the HJ model agree reasonably well with the experiments at all subsequent gauges, whereas using Mis-b005 yields similar results until 𝑥 = 12.35 m, but increasingly large differences in the crest area for shallower gauges. In contrast, the results of NWT-b005 yield surface elevations in better agreement with those predicted by Misthyc with the HJ energy dissipation model.

Periodic plunging breakers over a bar -Beji and

Battjes (1993) [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF], among others, performed laboratory experiments for periodic waves propagating over a trapezoidal bar (referred to here as BB-regular). Some of the generated incident waves were sufficiently steep to break over the bar as plunging breakers. Fig. 9 shows the set-up of the computational domain used in the models to simulate BB's experiments, with arrows indicating the locations of the wave gauge, where time series of the free surface elevation were measured in experiments. Waves were generated at 𝑥 = 0 in the constant depth region with ℎ 0 = 0.4 m, shoaled over the bar with a mild 1/20 offshore slope, and then broke over the crest of the bar where ℎ = 0.1 m, before the water depth increased again over the 1/10 onshore slope of the bar.

The target periodic wave characteristics in these experiments were a height 𝐻 0 = 0.054 m, period 𝑇 = 2.5 s, and wavelength 𝐿 0 = 4.8 m in the region of constant depth. However, in the wave gauge measurements, the actual wave height was 𝐻 0 ≃ 0.042 m at the toe of the slope, which was thus used as the wave height in the simulations. Simulations were performed with Misthyc and the NWT using 𝑏 = 0.05 and, in both models, breaking onset occurred at 𝑥 𝑏 ≃ 12.2 m, as compared to between 12 and 13 m in experiments. The breaking termination was specified at 𝐵 off = 0.3 in the models, which occurs before the water depth increases shoreward and causes wave deshoaling (end of cyan shaded region in Fig. 9). Fig. 10, compares the measured time series of free surface elevations to those computed in the models. They agree well in the shoaling region, for 6 < 𝑥 < 12 m, and in the breaking region, for 12 < 𝑥 < 14 m. Larger differences start occurring for 𝑥 > 14 m, but wave breaking on top of the bar and deshoaling for 14 < 𝑥 < 17 m are simulated well overall. Differences between the results of both models are small, except for higher harmonics generated on top of and beyond the bar, which Misthyc does not capture as well as the NWT. This may be due to the limited number of vertical layers (𝑁 𝑇 ), or the difference in timestep used.

The effect of the breaking termination criterion value was tested in the Misthyc model in the range 𝐵 off = 0.2 to 0.35. Fig. 11 shows the simulated and measured spatial variation of the significant wave height 𝐻 𝑠 throughout the domain. The results show that, in the zone after wave breaking (for 𝑥 > 14 m), using 𝐵 off = 0.20 leads to underpredicting 𝐻 𝑠 , whereas using 𝐵 off = 0.35 leads to overpredicting it. Using 𝐵 off = 0.30 or even 0.25, allows reproducing well 𝐻 𝑠 measured in the experiments after the bar. Thus, 𝐵 off = 0.30 was selected for all cases shown here involving submerged bars (where waves deshoal).

Irregular plunging breakers over a bar -Beji and Battjes (1993)

Using the same submerged bar, water depth, and set of wave gauges (see Fig. 9), [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF] also tested cases with irregular incident wave trains generated at the wavemaker based on a JONSWAP spectrum with significant wave height 𝐻 𝑠 = 0.049 m and peak spectral period 𝑇 𝑝 = 1∕𝑓 𝑝 = 2.5 s (referred to as BB-irregular). In the simulations, performed here using only Misthyc with 𝑏 = 0.05, incident waves are specified as a linear superposition of periodic waves in the frequency range [0.25𝑓 𝑝 , 5𝑓 𝑝 ], obtained from an FFT of the experimental free surface elevation measurements at the wave gauge located at 𝑥 = 6 𝑚. Once again, 𝐵 off = 0.30 is used in the model.

The simulation is run for 600 s or 240 peak periods, and breaking of the steepest waves was observed over the bar crest, as in the experiments where they were plunging breakers. To compare the simulated and experimental results, several wave statistics are computed based on the time series of free surface elevations at the wave gauges. These quantify the wave height (𝐻 𝑠 ) and wave shape and non-linear (1994) for regular spilling waves (circles) with simulations using Misthyc with the HJ model (dashed red line), Mis-b005 (solid red line), and NWT-b005 (solid blue line).

characteristics (asymmetry 𝐴𝑠, skewness 𝑆𝑘, kurtosis 𝐾𝑢), defined as follows,

𝐻 𝑠 = 4𝜎 1∕2 (18) 𝐴𝑠 = ⟨ℍ(𝜂 -⟨𝜂⟩) 3 ⟩∕𝜎 3∕2 (19) 𝑆𝑘 = ⟨(𝜂 -⟨𝜂⟩) 3 ⟩∕𝜎 3∕2 (20) 𝐾𝑢 = ⟨(𝜂 -⟨𝜂⟩) 4 ⟩∕𝜎 2 -3 (21)
where 𝜎 = ⟨(𝜂 -⟨𝜂⟩) 2 ⟩, is the free surface variance, ⟨ ⟩ the time averaging operator, and ℍ, the Hilbert transform.

Fig. 12 shows that there is a good agreement between the simulated and measured wave statistics for all wave characteristics. More specifically, Fig. 12a shows that 𝐻 𝑠 increases up to the onset of wave breaking at 𝑥 ≃ 12 m, then decreases during breaking over the bar crest (12 < 𝑥 < 14 m), and remains roughly constant for 𝑥 > 14 m. Fig. 12b shows that wave asymmetry is maximum at breaking onset (𝑥 ≃ 12 m) and decreases during wave breaking. Finally, Figs. 12c,d shows that the skewness and kurtosis, which quantify the wave nonlinearity, gradually increase during shoaling and breaking, and then decrease during deshoaling, after the bar. Both of these statistics agree better with the experiments than the results reported by [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], who also simulated this test case using Misthyc with the HJ dissipation and an eddy viscosity model [START_REF] Kurnia | High order hamiltonian water wave models with wave-breaking mechanism[END_REF] in combination with several breaking onset criteria that were not based on 𝐵. x (m)

-0.4 -0.2 0 z (m)
1/10 1/20

Figure 9: Set-up of the computational domains in the simulations of the [START_REF] Beji | Experimental investigation of wave propagation over a bar[END_REF] experiments for periodic waves shoaling and plunging breaking over a bar (BB). Arrows indicate the locations of wave gauges in the experiments and simulations. Wave breaking occurs in the models in the cyan shaded region, from onset to termination, and the gray shaded region indicates the absorbing beach. Note, the free surface shown is of BB-regular.

Irregular waves breaking over a bar -Adytia et al. (2018)

Adytia et al. ( 2018) also performed experiments with irregular waves propagating and breaking over a submerged bar with a different geometry than BB's (referred to as AH). Irregular incident waves were generated with 𝐻 𝑠 = 0.2 m and 𝑇 𝑝 = 2.5 s, in a constant depth ℎ 0 = 0.615 m. Fig. 13 shows the Misthyc model set-up, with the bar having a 1/20 offshore slope with the toe located at 𝑥 = 23.65 m and extending up to 𝑥 = 31.98 m, followed by a constant depth crest with ℎ = 0.2 m, to 𝑥 = 41 m. Time series of the free surface elevation were measured at 15 wave gauges, with the locations indicated with arrows in Fig. 13. As in the BBirregular test case, the free surface elevation measured at the gauge located at 𝑥 = 11.5 m is used to calculate the incident waves conditions for the simulation (based on a FFT).

The simulation is run for 500 s or 200 peak periods, and breaking onset occurs at 𝑥 𝑏 ∈ [30 -32] m. The same wave statistics as in the previous application are computed here based on experimental measurements and simulations with the Misthyc model, using 𝑏 = 0.05. These are shown in Fig. 14 where, overall, there is a good agreement between the experiments and simulations. Fig. 14a shows that 𝐻 𝑠 gradually increases over the offshore slope of the bar until wave breaking begins over the bar crest, then decreases in the breaking region (30 < 𝑥 < 35 m), and finally remains constant for 𝑥 > 35 m to the end of the bar crest. A similar trend is observed in asymmetry as in the BB-irregular test case in Fig. 14b: an increase up to the onset of breaking, a decrease in the breaking region, and a roughly constant value afterwards. Although the variation of skewness and kurtosis are reproduced better in the breaking region compared to [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], some differences are seen in the postbreaking region (see Figs. 14c,d).

Discussion

In the previous section, the results of simulations using one or both FNPF models with the newly proposed breaking onset/termination and dissipation parameterizations were presented, for standard experimental benchmarks from the literature (referred to as HS, BB, TK and AH) featuring regular or irregular waves propagating and breaking over several mild beach slopes and bars. In each case, the type of breaking (spilling-S or plunging-P) was reported in the experiments, as well as the measured breaking index value 𝜅 𝑏𝑒 = (𝐻 𝑏 ∕ℎ 𝑏 ) 𝑒 for some cases. Overall, the simulation results agreed well with the experimental data, confirming that the breaking onset criterion (𝐵 = 0.85) proposed by [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] is accurate, and wave elevation and kinematics at breaking onset and during breaking are accurately simulated in the models using the absorbing pressure and the proposed parameterizations of the energy dissipation rate.

Table 3 summarizes, for the periodic wave cases reported in the HS, BB-regular and TK studies, the incident wave and bathymetric parameters specified in the models and experiments: (𝐻 0 , 𝑇 ), offshore slope 𝑆, and wave characteristics simulated at breaking onset (𝑇 𝑏 , 𝑐 𝑏 , 𝜅 𝑏 ). Note the breaking period used in the definition of 𝛾, 𝑇 𝑏 < 𝑇 is based on the breaking crest geometry, following [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF]. In addition to the test cases described in detail in Section 4.2, additional simulations were performed for periodic wave cases from two more studies by [START_REF] Narayanan | Vertical structure of horizontal velocity in regular shoaling waves[END_REF] (NM) and [START_REF] Blenkinsopp | Void fraction measurements in breaking waves[END_REF] (BC), whose parameters and results are also summarized in Table 3. For each case, the Table also lists the value of Battjes' surf similarity parameter (or Irribaren number), 𝜉 0 = 𝑆∕ √ 𝐻 0 ∕𝐿 0 , where 𝐿 0 is the incident wavelength in deep water, 𝐿 0 = 𝑔𝑇 2 ∕(2𝜋) [START_REF] Battjes | Surf similarity[END_REF]. Battjes showed that periodic waves shoaling over a plane slope break as spilling breakers for 𝜉 0 ≤ 0.5 and as plunging breakers for 0.5 < 𝜉 0 ≤ 3.3. In all of the experiments reported in Table 3, waves broke either as spilling or plunging breakers, except in one case (c: S/P-BC) where both were observed. In most cases, the 𝜉 0 value is consistent with the observed type of wave breaking, despite some of these experiments being performed over a bar and not just a plane slope. Consistent with Battjes's work and predictions based on the surf similarity parameter for periodic depth-limited breakers, the models predicted a breaking index 𝜅 𝑏 ∈ [0.7 -1.2] for the tests considered here. In the next section, the dependence of the breaking strength 𝑏 predicted using the HJ model on the instantaneous 𝜅 ≃ 2𝐹 values in the breaking area is examined.

HJ model breaking strength

For the HJ model, the instantaneous breaking strength parameter 𝑏 can be computed in each case using Eq. ( 13) and expression of 𝜖 in Eq. ( 15), assuming Π 𝑏 = 𝜖. As shown in Eq. ( 14), assuming small amplitude waves (i.e., a small 𝐹 ) yields 𝑏 ≃ 2𝜇𝐹 3 with 𝐹 = 𝑔𝑎∕𝑐 2 𝓁 ≃ 𝑎∕𝑑 and 𝑐 𝓁 ≃ √ 𝑔𝑑 in shallow water. Since waves tend to not be symmetric about the mean water level, this equation underestimates the product ℎ 𝑐 ℎ 𝑡 in Eq. ( 13), and therefore overestimates the non-dimensional breaking strength compared to using the complete formula. Therefore, the value predicted this way represents an upper bound, i.e., 𝑏 𝑚𝑎𝑥 = 2𝜇𝐹 3 . Now, noting that in Eq. ( 13) we always have, ℎ 𝑐 ℎ 𝑡 < (𝑑 + 2𝑎) 2 = (1+2𝐹 ) 2 𝑑 2 , replacing the latter value in the equation yields a lower bound of the breaking strength, i.e., 𝑏 𝑚𝑖𝑛 = 𝑏 𝑚𝑎𝑥 ∕(1+2𝐹 ) 2 . As the energy dissipates during the breaking process, the relative depths under the crest and trough, ℎ 𝑐 and ℎ 𝑡 , change, and therefore the exact value of 𝑏 will vary, roughly bounded by these approximations, until the breaking termination criterion is met.

Fig. 15 shows the breaking strength 𝑏 computed with the Misthyc model, using the HJ parameterization, as a function of 𝜅(𝑡) = 𝐻∕𝑑 ≃ 2𝐹 , for the HS, TK, and BB-regular periodic wave experiments, as well as its approximate lower and upper bounds, 𝑏 𝑚𝑖𝑛 = (𝜇∕4)𝜅 3 ∕(1 + 𝜅) 2 and 𝑏 𝑚𝑎𝑥 = (𝜇∕4)𝜅 3 , respectively. Breaking onset is at the right side of the figure, where the largest values of 𝑏, 𝜅 and 𝐹 occur and, moving to the left of the figure, the waves propagate through the breaking region as a function of time, with breaking termination occurring near the bottom left. The 3).

These observations, as well as the good agreement between the simulations and experimental observations for a variety of benchmark cases, confirm the relevance of modeling the energy dissipation rate for depth-induced breaking waves as analogous to that of a hydraulic jump (or bore), with a single calibration constant 𝜇 = 1.5 [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF][START_REF] Stive | Energy dissipation in waves breaking on gentle slopes[END_REF] Eq. ( 13)). Furthermore, unlike in earlier studies that used a variety of case specific breaking criteria and related energy dissipation rates, good agreement with all experiments reported here for spilling and plunging breaker cases was achieved using universal breaking onset and termination criteria based on 𝐵 and a dissipation rate based on the HJ analogy. Fig. 15 also shows that the timeaveraged breaking strength, based on all simulated cases, is 𝑏 ≃ 0.05, as tested in several simulations in Section 4.2. And although the actual 𝑏 value differs substantially from this average near breaking onset and termination, as seen in earlier results, there were no significant differences in the results obtained with a constant 𝑏 = 0.05 or varying 𝑏 value throughout breaking within the accuracy/variability of most experiments, and both approaches agreed similarly with the experiments for the 5 considered benchmark cases.

Breaking onset kinematics

Fig. 16 shows the time evolution of the 𝐵 value up to breaking onset (𝐵 = 0.85), as computed in simulations of the 10 experimental cases listed in Table 3. As expected from the parameterization of 𝑏 as an increasing function of 𝛾 ∝ 𝑑𝐵∕𝑑𝑡, defined in Eq. ( 16) (Fig. 3), 𝑑𝐵∕𝑑𝑡 is observed to be larger at breaking onset for plunging breakers than for the spilling breakers. This observation appears to be independent from the type of bathymetry that caused waves to break. For instance, BB and TK conducted S and P breaker experiments on barred and sloped bathymetries, respectively, and 𝑑𝐵∕𝑑𝑡 is larger for P than for S breakers in all cases. [START_REF] Blenkinsopp | Void fraction measurements in breaking waves[END_REF] types of experiments, for strongly P (BC1), P (BC2), and S/P (BC) breakers and, among these, 𝑑𝐵∕𝑑𝑡 is largest for BC1 and lowest for BC, as would be expected.

For the cases shown in Table 3, the S breakers, with 𝜉 0 < 0.7 have 𝛾 < 1.3 and the P breakers, have 𝜉 0 > 0.5 and 𝛾 > 1.4. Hence, consistent with the value of 𝑑𝐵∕𝑑𝑡 at breaking onset, 𝛾 values may also distinguish S from P breakers. However, a more accurate estimation of this threshold 𝛾 value should be based on a larger number of test cases.

Breaking termination conditions

For waves breaking on a plane beach (such that depth always decreases as waves propagate into shallower water), as described in [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF], from breaking onset shoreward, the breaking dissipation can occur up until waves run up the dry upper slope. However, the FNPF setup used here does not model wave run-up or bottom friction, which becomes significant near the shoreline. Using a similar model set-up, simulations in previous work on plane beaches did not use a breaking termination criterion. For example, in Fig. 9a from [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], the wave height decreases from the onset up to the shallowest depth. For the cases on a plane slope, this is also seen in Fig. 10 of [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], where the same variation in 𝐻 𝑠 is observed. Finally, [START_REF] Papoutsellis | Modelling of depth-induced wave breaking in a fully nonlinear free-surface potential flow model[END_REF] also mention that for the spilling breaker case of [START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF], they terminate breaking inside the sponge layer.

In contrast, for waves breaking over a bar, which is followed by deshoaling over the shoreward slope of the bar, a breaking termination criterion 𝐵 off = 0.30 was used in all of the present applications. As shown in Fig. 11, this value appears to be optimal for reproducing the experimental results. Waves breaking on a bar may terminate breaking over the crest or on the shoreward slope of the bar. Here, the choice of termination criterion is important since the wave evolution depends strongly on it.

To the authors' knowledge, a universal breaking termination criterion does not exist. Alternative ideas have been tested, such as a termination criterion based on the slope of the free-surface (e.g., [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF]. Another possibility is a time-based criterion. Derakhti et al. (2018a) noted that the active breaking period 𝜏 ≈ 0.75𝑇 𝑏 , and so this could also be tested in future applications.

Conclusions

A unified method of modeling depth-limited wave breaking dissipation in FNPF models was demonstrated, building on the work of [START_REF] Guignard | Modeling of wave shoaling in a 2D-NWT using a spilling breaker model[END_REF] and [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF], who proposed making the energy dissipation rate analogous to that of a hydraulic jump. Two different FNPF models were applied here: Misthyc [START_REF] Yates | Accuracy and efficiency of two numerical methods of solving the potential flow problem for highly nonlinear and dispersive water waves[END_REF] and a BEM-NWT [START_REF] Grilli | An efficient boundary element method for nonlinear water waves[END_REF][START_REF] Grilli | Numerical modeling of wave breaking induced by fixed or moving boundaries[END_REF]. The detection of breaking onset in the models was based on the universal criterion 𝐵 = 𝑢∕𝑐 = 0.85, first proposed by [START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF] for deep and intermediate water cases, and validated by [START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF] in shallow water. Similar to earlier work such as [START_REF] Grilli | Fully nonlinear potential flow simulations of wave shoaling over slopes: Spilling breaker model and integral wave properties[END_REF] or [START_REF] Simon | Comparing methods of modeling depth-induced breaking of irregular waves with a fully nonlinear potential flow approach[END_REF], we show that simulation results based on this approach agree well with experimental measurements for a variety of standard shallow water breaking cases from the literature, for both regular and irregular wave trains and different bathymetries.

Using the hydraulic jump (HJ) analogy originally proposed by [START_REF] Svendsen | Wave characteristics in the surf zone[END_REF] and [START_REF] Stive | Energy dissipation in waves breaking on gentle slopes[END_REF] for spilling breakers, the non-dimensional breaking strength 𝑏, is found in most cases, including both spilling or plunging breakers, to have a fairly narrow range of variation centered on 0.05. However, consistent with the recent parameterization of 𝑏 proposed by Derakhti et al. (2018a,b), proportional to 𝑑𝐵∕𝑑𝑡 at the onset, the instantaneous value of 𝑏 resulting from the HJ analogy is much larger at breaking onset, and then decreases throughout the breaking region. Recalling that the HJ dissipation rate is related to the relative wave height 𝜅, or wave Froude number 𝐹 , with 𝜅 = 𝐻∕𝑑 ≃ 2𝐹 in shallow water and for the depth-limited breaking waves tested here, 𝜅 𝑏 ∈ [0.7 1.2], consistent with Battjes' predictions based on the surf similarity parameter [START_REF] Battjes | Surf similarity[END_REF]. This implies that for waves in this parameter range, 𝑏 would always have a similar order of magnitude. With this rationale, simulations of spilling and plunging depth-limited breaking waves were performed using a constant breaking strength 𝑏 = 0.05. The simulated results using a constant 𝑏 or the HJ model with a time-varying 𝑏 agreed similarly with the experimental measurements, within the range of experimental uncertainty. Although in very complex cases there may be waves in an irregular wave train requiring larger or smaller 𝑏 values, the constant 𝑏 value provides results with a similar level of uncertainty or accuracy as the HJ model.

There are some weaknesses in this current approach that must be resolved in order to obtain a fully general method. While concentrating on the wave breaking criterion and dissipation rate, the breaking termination criterion has not yet been thoroughly investigated, which should be a point for further research. For the cases shown here, wave breaking on a slope were not terminated until the shallowest depth in the model domain (or the absorbing beach) were reached, whereas wave breaking on a bar was terminated at a value 𝐵 off = 0.30, which was found to yield results in good agreement with experiments. Since the energy dissipation rate Π 𝑏 = 𝜖 of breaking waves depends on the fifth power of wave celerity, when this parameterization is used in the constant strength model, an accurate tracking of wave crest locations as a function of time 𝑥 𝑐 (𝑡) and the calculation of the time derivative 𝑑𝑥 𝑐 ∕𝑑𝑡 = 𝑐 are important for accurate simulations. With the method used here to compute 𝑐 (see Appendix), some spurious oscillations are still observed when taking the derivative of the crest position, which results in small oscillations of Π 𝑏 = 𝜖 and thus wave heights in the breaking region. Another issue is the accurate realtime identification of wave crests/troughs in irregular seastates. In situations where small breaking waves ride on top of larger waves, the estimation of the spatial extent of dissipation could be misinterpreted, which could cause instabilities.

By using a breaking dissipation rate Π 𝑏 = 𝜖 proportional to a constant 𝑏, instead of using the HJ analogy, the need to compute the local geometric characteristics of waves required in the HJ dissipation Eq. ( 13) is eliminated (e.g., water depth under the trough, etc.), and the dissipation can simply be related to the crest kinematics. This simpler formulation clearly will yield its greatest advantages when applying the model in three-dimensions, where the identification of local wave parameters, even the crest location, is challenging. Derakhti et al. (2018b[START_REF] Derakhti | Predicting the rate of energy dissipation in breaking gravity water waves from deep to shallow water[END_REF]) also showed that a constant 𝑏 might be expected in shallow water based on results of a Navier-Stokes model. The present modeling approach with the constant strength 𝑏 and dissipation rate Π 𝑏 = 𝜖 parameterization could be generalized to intermediate or deep water cases, by using the 𝑏(𝛾) breaking strength model proposed by Derakhti et al. (2018a), who found that the value of 𝑏 in deep or intermediate water (which is clearly not constant) could be related to the kinematics of the wave as well, specifically to 𝑑𝐵∕𝑑𝑡 and a relevant breaking period 𝑇 𝑏 at breaking onset. Some additional work may be required, however, to implement the deep-water parameterization, as 𝑏 obtained from a 3D finite volume NS solver, in which energy is dissipated within the domain volume, may not have the same distribution or effect in a FNPF model, where the energy dissipation is applied to the surface and only acts on the potential part of the flow.

In considering the kinematics of the shallow water breaking waves simulated here, the value of 𝛾 appears to provide a parameter to distinguish spilling and plunging breakers, with a threshold value around 1.3 -1.4. For waves shoaling on a plane slope, this may be loosely related to Battjes' surf similarity parameter and may be of interest for applications beyond the scope of those envisioned here, such as of identifying breaker types from existing models that are unable to simulate the breaking process.

A. Calculation of wave crest kinematics

The instantaneous location of wave crests is first roughly estimated on the 2D-FNPF models' free surface, by finding local maxima. Then, four model points surrounding each wave crest location are identified with elevation, particle velocity, and location (𝜂 𝑖 , 𝑢 𝑖 , 𝑥 𝑖 ) (𝑖 = 1, 2, 3, 4) and mapped to a local coordinate 𝜉 ∈ [-1, 1]. Cubic shape functions (Eq. 22) are used to interpolate 𝜂(𝜉) between these points.

The wave crest location is then estimated based on the 𝜉 𝑐 value at which 𝑑𝜂∕𝑑𝜉 = 0 (Fig. 17). Therefore, the horizontal location of a local crest in the global coordinate is calculated as 𝑥 𝑐 = 𝑥 1 𝑁 1 (𝜉 𝑐 ) + 𝑥 2 𝑁 2 (𝜉 𝑐 ) + 𝑥 3 𝑁 3 (𝜉 𝑐 ) + 𝑥 4 𝑁 4 (𝜉 𝑐 ), (23) and the horizontal particle velocity at the crest, 𝑢 as 𝑢 = 𝑢 1 𝑁 1 (𝜉 𝑐 ) + 𝑢 2 𝑁 2 (𝜉 𝑐 ) + 𝑢 3 𝑁 3 (𝜉 𝑐 ) + 𝑢 4 𝑁 4 (𝜉 𝑐 ). where (𝛼, 𝛽) based on some preliminary tests are taken as (0.05, 0.01).

Wave celerity of a crest at a time 𝑡, 𝑐(𝑡) is then calculated by a linear fit to the smoothed data in the interval [𝑋 𝑡-𝑛 𝑐 , ..., 𝑋 𝑡 𝑐 ], such that 𝑛Δ𝑡 < 0.04𝑇 , where Δ𝑡 is the discretization in time and 𝑇 is a representative wave period. In the test cases presented this study, 𝑛 = 8. To compute 𝛾, the discrete values of 𝐵 = 𝑢∕𝑐 at the crest, when in the range [0.82 -0.85], are then linearly fit to obtain 𝑑𝐵∕𝑑𝑡| 𝐵=𝐵 𝑡ℎ .

Figure 1 :

 1 Figure 1: (top) Geometric parameters used to calculate the HJ dissipation, and (bottom) the shape function 𝑆(𝑥)

Figure 2 :

 2 Figure 2: Evolution of 𝐵 = 𝑢∕𝑐 in simulations with the Misthyc model as a function of non-dimensional time 𝑡 * = (𝑡 -𝑡 𝑏 )∕𝑇 𝑏 (𝑡 𝑏 denotes the time of breaking onset when 𝐵 = 𝐵 𝑡ℎ = 0.85), for periodic: (i) plunging breaking waves propagating over a bar, from Beji and Battjes (1993) (BB-regular, red); (ii) spilling breaking waves propagating over a plane slope, from Hansen and Svendsen (1979) (HS, magenta); and shoaling and spilling breaking waves propagating over a plane slope, from Ting and Kirby (1994) (TK, green).

  (BB-regular); (ii) shoaling and spilling breaking waves propagating over a plane slope, from Hansen and Svendsen (1979) (HS); and (iii) shoaling and spilling breaking waves propagating over a plane slope from Ting and Kirby (1994) (TK). The details of the set-up and numerical parameters for these simulations are described in the following subsections. Note, the HJ dissipation model and corresponding Π 𝑏 values were previously experimentally validated using the NWT model for the HS test cases by Grilli et al. (2020), and using the Misthyc model for the regular and irregular BB cases by Simon et al. (2019), and with another FNPF model called HCMT, for the regular wave TK and BB cases by Papoutsellis et al. (2019).

Fig. 2

 2 shows the evolution of 𝐵 computed leading up to wave breaking with the Misthyc model as a function of the non-dimensional time 𝑡 * = (𝑡 -𝑡 𝑏 )∕𝑇 𝑏 , where 𝑡 𝑏 is the time of breaking onset over the region of interest for calculating 𝑑𝐵∕𝑑𝑡. Consistent withDerakhti et al. (2018a,b), of the three test cases shown, the plunging breaker case (BBregular) has the fastest rate of change of 𝐵 at the threshold, whereas the spilling breaker cases (HS and TK) show slower changes in 𝐵. The calculation of 𝑑𝐵∕𝑑𝑡 was found to be sensitive to the interval over which it is calculated, as well as to the spatial and temporal resolution of the simulation. Fig.2(right) shows the linear fit applied to 𝐵 in the interval [0.82-0.85], used to compute 𝑑𝐵∕𝑑𝑡 and, for this interval for the HS case, Table1shows the sensitivity of the 𝑑𝐵∕𝑑𝑡 and 𝛾 values calculated for different spatial and temporal grids. Fig. 3 shows the average 𝑏 values computed in each test case using the HJ model as a function of 𝛾, compared to the values predicted by the second parameterization based on 𝛾 (Eq. 16). Average values of 𝑏 computed over the breaking region using the HJ dissipation range within [0.02, 0.08], and corresponding 𝛾 vary within [0.6, 2.5]. These results are in moderate agreement with those of the Derakhti et al. (2018a,b) parameterization, when 𝛾 ∈ [1 -2]. Some differences between these two parameterizations of 𝑏 are to be expected since even the Derakhti et al. (2018a,b) parameterization is only a curve fit, with significant spread of individual values.

Figure 3 :

 3 Figure 3: Average breaking strength 𝑏 computed for the HJ parameterization, as a function of 𝛾 for depth-limited breaking waves simulated in Misthyc including: (filled circles) regular wave cases, magenta: HS, green: TK and red: BB; (hollow circles) BB-irregular cases. For comparison, the empirical parmaterization from Eq. 16 proposed by Derakhti et al. (2018a), is indicated with a solid line. The dashed line shows the constant strength average breaking value 𝑏 = 0.05 used in the Mis-005 parameterization.

Figure 4 :

 4 Figure 4: Comparison of the instantaneous energy dissipation rate Π 𝑏 (𝑡) calculated for the plunging BB-regular case in the Misthyc model using: (blue) the HJ analogy Eq. (13), or (red) the Mis-b005 Eq. (15) parameterization with 𝑏 = 0.05.

Figure 5 :

 5 Figure 5: Sensitivity of the simulation results to the spatial discretization for the HS test case using the NWT (blue dashed, 𝐿∕Δ𝑥 = 35, and blue solid, 𝐿∕Δ𝑥 = 50) and Misthyc (solid red with 𝐿∕Δ𝑥 = 70) model.

Figure 6 :

 6 Figure 6: (a) Bathymetry of the Hansen and Svendsen (1979) (HS) experimental set-up for periodic shoaling and spilling breaking waves, where the shaded cyan region indicates the wave crest location from the onset of wave breaking (note an absorbing beach is specified in the models for 𝑥 ≥ 25.5 m). Spatial evolution of the (b) wave celerity, and (c) wave height in the experimental data (circles), Misthyc using the HJ model (dashed red line), Mis-b005 (solid red line) and NWT-b005 (solid blue line), averaged over 5 successive wave periods after a quasi-steady state is reached.

Figure 7 :

 7 Figure 7: Setup for Ting and Kirby (1994) laboratory experiments of periodic spilling breakers propagating over a 1/35 slope, with black arrows showing the locations of the 12 wave gauges. The cyan shading indicates the simulated wave breaking region, from onset to termination. The gray shading indicates the beach absorption zone.

Figure 8 :

 8 Figure8: Comparison of the temporal evolution of the free surface position measured at the gauge locations from[START_REF] Ting | Observation of undertow and turbulence in a laboratory surf zone[END_REF] for regular spilling waves (circles) with simulations using Misthyc with the HJ model (dashed red line), Mis-b005 (solid red line), and NWT-b005 (solid blue line).

Figure 10 :

 10 Figure 10: Comparison of the temporal evolution of the free surface at the gauge locations from the Beji and Battjes (1993) test case for regular plunging waves (circles) with Mis-b005 (red solid line) and NWT-b005 (blue solid line).

Figure 11 :

 11 Figure 11: Spatial evolution of the significant wave height 𝐻 𝑠 for the BB-regular experiments (circles) and simulations with the Mis-b005 model using different breaking termination criteria 𝐵 off =: (red) 0.2, (blue) 0.25, (green) 0.30, and (yellow) 0.35.

Figure 12 :

 12 Figure 12: Spatial evolution of wave statistics computed from the experiments (circles) and simulations with Mis-b005 (solid line) for the irregular wave case of Beji and Battjes (1993) (JONSWAP spectrum with 𝐻 𝑠 = 0.049 m and 𝑇 𝑝 = 2.5 s). Waves break over the bar in the model and in the experiments (plunging breakers).

Figure 13 :

 13 Figure 13: Set-up for simulations with Misthyc of Adytia et al. (2018) experiments for irregular waves propagating over a bar (𝐻 𝑠 = 0.2 m and 𝑇 𝑝 = 2.5 s), with arrows showing the location of the 15 wave gauges. The cyan shading shows the breaking region from onset to termination, and the gray shading indicates the absorbing beach region.

Figure 14 :

 14 Figure 14: Spatial evolution of wave statistics computed based on results of experiments (circles) and simulations with Mis-b005 (solid line) for the irregular wave case of Adytia et al. (2018) (𝐻 𝑠 = 0.2 m and 𝑇 𝑝 = 2.5 s); waves break over the bar in the model and in the experiments.

Figure 15 :

 15 Figure 15: Evolution of the instantaneous breaking strength 𝑏 calculated with the HJ model during breaking simulations with Misthyc for the (Table3): (a) HS (solid magenta), (b) TK (solid green), and (i) BB-regular (solid red) periodic spilling/plunging breaking wave experiments, as a function of 𝜅 = 𝐻∕𝑑. The dashed blue curves represent the lower (𝑏 𝑚𝑖𝑛 ) and upper (𝑏 𝑚𝑎𝑥 ) bounds estimated for 𝑏, and the horizontal dash-dotted black line is 𝑏 = 0.05.

Figure 16 :

 16 Figure 16: Evolution of 𝐵 = 𝑢∕𝑐 as a function of non-dimensional time 𝑡 * = (𝑡 -𝑡 𝑏 )∕𝑇 𝑏 , referred to the time of breaking onset, up to breaking onset (𝐵 = 0.85; horizontal dashed line), computed for a wave crest evolving over sloping bathymetries in the experimental test cases listed in Table3. Red lines for P-cases and blue lines for S-cases have smaller and larger values of 𝑑𝐵∕𝑑𝑡 at breaking onset, respectively.

Figure 17 :

 17 Figure 17: Sketch showing the calculation of a local wave crest location between discretization points on the simulated free surface.

  (24) Prior to calculating 𝑐, 𝑥 𝑐 (𝑡) is smoothed by second-order exponential smoothing[START_REF] Guthrie | NIST/SEMATECH e-Handbook of Statistical Methods[END_REF], described as follows, denoting 𝑋 𝑡 𝑐 as the smoothed horizontal crest location at time𝑡 = 0, 𝑋 0 𝑐 = 𝑥 0 𝑐 𝑠 0 = 𝑥 1 𝑐 -𝑥 0 𝑐 (25)For 𝑡 > 0,𝑋 𝑡 𝑐 = 𝛼𝑥 𝑡 𝑐 + (1 -𝛼)(𝑋 𝑡-1 𝑐 + 𝑠 𝑡-1 ) 𝑠 𝑡 = 𝛽(𝑋 𝑡 𝑐 -𝑋 𝑡-1 𝑐 ) + (1 -𝛽)𝑠 𝑡-1 ,

  [START_REF] Derakhti | Wave breaking in the surf zone and deep-water in a non-hydrostatic RANS model. part 1: Organized wave motions[END_REF] extended Duncan's parameterization of the rate of energy dissipation 𝜖 to intermediate water, steepness-limited, plunging breaking waves and irregular wave trains. Following the introduction of the breaking onset criterion based on 𝐵 by[START_REF] Barthelemy | On a unified breaking onset threshold for gravity waves in deep and intermediate depth water[END_REF],Derakhti et al. (2018a) observed in their NS simulations of focused wave trains, that the slope of 𝐵(𝑡) at breaking onset was correlated with the average breaking power dissipated in the model, from breaking onset to termination, ⟨𝜖⟩, or the corresponding breaking strength 𝑏 obtained from Eq. (15). On this basis, they proposed a new parameterization of the breaking strength parameter, 𝑐 𝑙𝑏 = √ 𝑔 tanh(𝑘 𝑏 𝑑)∕𝑘 𝑏 , with 𝑘 𝑏 = 2𝜋∕𝐿 𝑏 , the linear wave celerity at breaking in arbitrary depth, and 𝐿 𝑏 a relevant breaking wave length calculated based on the method proposed by

	𝑏 =	𝑔⟨𝜖⟩ 𝑙𝑏 𝜌𝑐 5	= 0.034 (𝛾 -0.30) 2.5	(16)
		𝛾 = 𝑇 𝑏	𝑑𝐵 𝑑𝑡	| |𝐵=𝐵 𝑡ℎ | |	(17)

computed based on wave crest parameters at breaking onset and a breaking wave period, 𝑇 𝑏 = 𝐿 𝑏 ∕𝑐 𝑙𝑏 , estimated based on,

Table 1 Sensitivity

 1 Δ𝑥 (𝑚) Δ𝑡 (𝑠) 𝑑𝐵𝑑𝑡 | 𝐵=𝐵 𝑡ℎ (𝑠 -1 ) 𝛾

	0.028	0.007	0.885	0.726
	0.020	0.008	0.974	0.788
	0.031	0.008	1.00	0.774
	0.041	0.010	1.075	0.791

of 𝛾 computed with Eq. (

16

) for different spatial and temporal resolutions, for the

[START_REF] Hansen | Regular Waves in Shoaling Water Experimental Data[END_REF]

(HS) case.

Table 2

 2 Numerical parameters used in the simulations with the Misthyc and NWT models for the five experimental benchmark test cases.

	Test case	Δ𝑥 (𝑚)		Δ𝑡 (𝑠)	Domain	𝑇 𝑚𝑎𝑥
		Misthyc	NWT	Misthyc	NWT	length (m)	(𝑠)
	HS	0.019	0.031	0.008	0.008	22.7	25
	TK	0.019	0.055	0.008	0.013	29.7	32
	BB-regular	0.014	0.048	0.01	0.015	34.4	35
	BB-irregular	0.02	-	0.01	-	30	600
	AH	0.02	-	0.01	-	60	500

  conducted three

	Label	Name	𝐻 0 (cm) 𝑇 (s) slope 𝑆 𝑇 𝑏 (s) 𝑐 𝑏 (m/s)	𝜅 𝑏	𝜅 𝑏𝑒	𝛾	𝜉 0
	a	M-S-HS	9.50	1.0	1/34.26	0.78	1.03	0.73 0.78 0.77 0.118
	b	M-S-TK	12.50	2.0	1/35	1.33	1.26	0.81 0.82 0.80 0.202
	c	N-S/P-BC	10.05	1.0	1/10	0.75	1.06	0.69 0.87 1.01 0.394
	d	N-S-NM	12.00	2.0	1/34.66	1.23	1.31	0.84	-	1.12 0.208
	e	N-P-NM	12.00	5.0	1/34.66	1.54	1.24	1.09	-	1.54 0.520
	f	M-S-BB	4.40	2.5	1/20	0.82	0.89	0.79	-	1.30 0.744
	g	N-P-TK	12.80	5.0	1/35	1.71	1.27	1.09 1.21 1.65 0.499
	h	N-P-BC2	9.72	1.42	1/10	1.04	1.07	0.95 0.91 1.43 0.569
	i	M-P-BB	5.40	2.5	1/20	0.99	0.92	0.92	-	2.47 0.672
	j	N-P-BC1	7.82	2.0	1/10	1.64	1.04	1.17 0.84 1.52 0.894
	Table 3
	Periodic wave breaking test cases and their parameters. Each test case name has 3 parts: (1) the model used, M-Misthyc, N-NWT;
	(2) type of breaking reported in experiments, S-Spilling, P-Plunging; and (3) the experiment: (slope) HS-Hansen and Svendsen
	(1979), TK-Ting and Kirby (1994), NM-Narayanan and McCalpin (1997) and BC-Blenkinsopp and Chaplin (2007), and (bar)
	BB-Beji and Battjes (1993).

Acknowledgements

This research was produced within the framework of Energy4Climate Interdisciplinary Center (E4C) of IP Paris and Ecole des Ponts ParisTech. This research was supported by 3rd Programme d'Investissements d'Avenir [ANR-18-EUR-0006-02]. This action benefited from the support of the Chair « Challenging Technology for Responsible Energy » led by l'X -Ecole Polytechnique and the Fondation de l'Ecole Polytechnique, sponsored by TOTAL.

SG is gratefully acknowledging support from the US National Science Foundation grant #OCE-19-47960.

The authors also thank Luc Pastur (ENSTA Paris) and Christophe Peyrard (EDF R&D LNHE) for helpful discussions.