Implementation and Evaluation of Breaking Detection Criteria for a Hybrid Boussinesq Model

TL;DR

This paper develops a hybrid wave model combining non-linear shallow water and Boussinesq equations, introducing a new wave breaking detection criterion based on free surface Froude number, with a novel discretization approach and extensive benchmarking.

Contribution

It presents a new physically meaningful wave breaking detection method and integrates it into a hybrid Boussinesq model with a novel finite element discretization strategy.

Findings

The new detection criterion effectively identifies wave breaking events.

The hybrid model accurately simulates wave propagation, shoaling, and run-up.

Benchmark tests demonstrate the robustness and advantages of the proposed method.

Abstract

The aim of the present work is to develop a model able to represent the propagation and transformation of waves in nearshore areas. The focus is on the phenomena of wave breaking, shoaling and run-up. These different phenomena are represented through a hybrid approach obtained by the coupling of non-linear Shallow Water equations with the extended Boussinesq equations of Madsen and Sorensen. The novelty is the switch tool between the two modelling equations: a critical free surface Froude criterion. This is based on a physically meaningful new approach to detect wave breaking, which corresponds to the steepening of the wave's crest which turns into a roller. To allow for an appropriate discretization of both types of equations, we consider a finite element Upwind Petrov Galerkin method with a novel limiting strategy, that guarantees the preservation of smooth waves as well as the monotonicity of the results in presence of discontinuities. We provide a detailed discussion of the implementation of the newly proposed detection method, as well as of two other well known criteria which are used for comparison. An extensive benchmarking on several problems involving different wave phenomena and breaking conditions allows to show the robustness of the numerical method proposed, as well as to assess the advantages and limitations of the different detection methods.

Equations

$$\begin{array}{rl}& \{\begin{array}{l}{\partial }_t\eta +{\partial }_{x}q=0\\ {\partial }_{t}q-\mathcal{D}(\eta ,q)+{\partial }_x(uq+g{H}^2/2)+gH{\partial }_{x}h+C_{f}q=0\end{array}\\ & \text{with}\mathcal{D}(\eta ,q)=B{h}^2{\partial }_{xxt}q+\beta g{h}^3{\partial }_{xxx}\eta +h{\partial }_{x}h{\partial }_{xt}q/3+2\beta g{h}^2{\partial }_{x}h{\partial }_{xx}\eta \end{array}$$

$$\{\begin{array}{l}{\partial }_t\eta +{\partial }_{x}q=0\\ {\partial }_{t}q-f_{\text{break}}(x,t)\mathcal{D}(\eta ,q)+{\partial }_x(uq+g{H}^2/2)+gH{\partial }_{x}h+C_{f}q=0,\end{array}$$

$$F(\mathbf{W})=\left[\begin{array}[]{c}q\\ uq+g\dfrac{H^{2}}{2}\end{array}\right]\!,\;S_{h}=gH\left[\begin{array}[]{c}0\\ \partial_{x}h\end{array}\right]\!,\;S_{f}=\left[\begin{array}[]{c}0\\ C_{f}q\end{array}\right]\!,\;D=\left[\begin{array}[]{c}0\\ \mathcal{D}(\mathbf{W})\end{array}\right].$$

$${\mathbf{W}}_{\text{h}}(t,x)=\sum _{j=0}^N{\phi }_j(x){\mathbf{W}}_j(t),F_{\text{h}}(t,x)=\sum _{j=0}^N{\phi }_j(x)F({\mathbf{W}}_j(t)).$$

$${\text{r}}_{\text{h}}={\partial }_t{\mathbf{W}}_{\text{h}}-f_{\text{break}}D({\mathbf{W}}_{\text{h}})+{\partial }_x{F}_{\text{h}}+S_h({\mathbf{W}}_{\text{h}})+S_f({\mathbf{W}}_{\text{h}}),$$

$$\begin{array}{r}\underset{\Omega }{\int }{\phi }_i{\text{r}}_{\text{h}}dx+\sum _{i=0}^{N-1}\underset{x_i}{\overset{x_i+\Delta x}{\int }}\left(A({\mathbf{W}}_{\text{h}}){\partial }_x{\phi }_i{\tau }_{\text{s}}\right){\text{r}}_{\text{h}}dx=0,\end{array}$$

$$\begin{split}\sum\limits_{l=-1,0,1}&M_{i,i+l}\Big{(}\dfrac{d\mathbf{W}_{i+l}}{dt}+S_{f_{i+l}}\Big{)}+\dfrac{\mathrm{I}-\mathrm{sign}\big{(}A(\mathbf{W}_{i+1/2})\big{)}}{2}\phi_{i+1/2}\\ &\qquad\qquad\qquad\qquad\qquad\;\;\;+\dfrac{\mathrm{I}+\mathrm{sign}\big{(}A(\mathbf{W}_{i-1/2})\big{)}}{2}\phi_{i-1/2}=\\ f_{\text{break}_{i}}&\int\limits_{\Omega}\varphi_{i}D(\mathbf{W}_{\textrm{h}})dx+\sum_{i=0}^{N-1}f_{\text{break}_{i+1/2}}\int\limits_{x_{i}}^{x_{i+1}}\left(A(\mathbf{W}_{\textrm{h}})\partial_{x}\varphi_{i}\,\tau_{\textrm{s}}\right)D(\mathbf{W}_{\textrm{h}})dx,\end{split}$$

$$\begin{split}\int\limits_{\Omega}\varphi_{i}\mathcal{D}(\mathbf{W}_{\textrm{h}})dx=\int\limits_{x_{i-1}}^{x_{i+1}}\big{(}-&\partial_{x}(Bh^{2}_{\textrm{h}}\varphi_{i})\partial_{xt}q_{\textrm{h}}+\varphi_{i}\dfrac{1}{3}h_{\textrm{h}}\partial_{x}h_{\textrm{h}}\partial_{xt}q_{\textrm{h}}\\ +&\beta g\varphi_{i}(h_{\textrm{h}}^{3}\partial_{x}w^{\eta}_{\textrm{h}}+2h_{\textrm{h}}^{2}\partial_{x}h_{\textrm{h}}\,w^{\eta}_{\textrm{h}})\big{)}dx\,,\end{split}$$

$$\begin{split}\int\limits_{x_{i}}^{x_{i+1}}&\left(A(\mathbf{W}_{\textrm{h}})\partial_{x}\varphi_{i}\,\tau_{\textrm{s}}\right)D(\mathbf{W}_{\textrm{h}})dx=-\frac{\text{sign}\big{(}A(\mathbf{W}_{i+\frac{1}{2}})\big{)}}{2}D_{i+1/2},\quad\text{with}\\ D_{i+1/2}=\int\limits_{x_{i}}^{x_{i+1}}&\left[\begin{array}[]{c}0\\[10.0pt] Bh_{\textrm{h}}^{2}\partial_{xt}w^{q}_{\textrm{h}}+\beta gh_{\textrm{h}}^{3}\partial_{x}w^{\eta}_{\textrm{h}}+h_{\textrm{h}}\partial_{x}h_{\textrm{h}}\partial_{xt}q_{\textrm{h}}/3+2\beta gh^{2}_{\textrm{h}}w^{\eta}_{\textrm{h}}\end{array}\right]dx.\end{split}$$

$$\underset{\Omega }{\int }{\phi }_i{w}_{\text{h}}^{\eta }=-\underset{i-1}{\overset{i+1}{\int }}{\partial }_x{\eta }_{\text{h}}{\partial }_x{\phi }_{i}dx,\underset{\Omega }{\int }{\phi }_i{w}_{\text{h}}^q=-\underset{i-1}{\overset{i+1}{\int }}{q}_{\text{h}}{\partial }_x{\phi }_{i}dx.$$

$$\begin{split}\Delta x\frac{d\mathbf{W}_{i}}{dt}+\dfrac{\mathrm{I}-\mathrm{sign}\big{(}A_{i+1/2}\big{)}}{2}\phi_{i+1/2}+\dfrac{\mathrm{I}+\mathrm{sign}\big{(}A_{i-1/2}\big{)}}{2}\phi_{i-1/2}=0\end{split}$$

$$\begin{split}\Delta x\frac{d\mathbf{W}_{i}}{dt}&+\frac{\Delta x\,\delta_{i-\frac{1}{2}}}{2}\Big{\{}\frac{1}{3}[\frac{d\mathbf{W}_{i-1}}{dt}-\frac{d\mathbf{W}_{i}}{dt}]+\text{sign}(A_{i-\frac{1}{2}})\dfrac{d\mathbf{W}_{i-\frac{1}{2}}}{dt}\Big{\}}\\ &+\frac{\Delta x\,\delta_{i+\frac{1}{2}}}{2}\Big{\{}\frac{1}{3}[\frac{d\mathbf{W}_{i+1}}{dt}-\frac{d\mathbf{W}_{i}}{dt}]-\text{sign}(A_{i+\frac{1}{2}})\dfrac{d\mathbf{W}_{i+\frac{1}{2}}}{dt}\Big{\}}\\ &+\dfrac{\mathrm{I}-\mathrm{sign}\big{(}A_{i+1/2}\big{)}}{2}\phi_{i+1/2}+\dfrac{\mathrm{I}+\mathrm{sign}\big{(}A_{i-1/2}\big{)}}{2}\phi_{i-1/2}=0\end{split}$$

$${\left({\partial }_{xx}\eta \right)}^{\text{O4}}=-\frac{1}{12\Delta x^2}{\eta }_{i+2}+\frac{4}{3\Delta x^2}{\eta }_{i+1}-\frac{5}{2\Delta x^2}{\eta }_i+\frac{4}{3\Delta x^2}{\eta }_{i-1}-\frac{1}{12\Delta x^2}{\eta }_{i-2}$$

$${\left({\partial }_{xx}\eta \right)}^{\text{O2}}=\frac{1}{\Delta x^2}{\eta }_{i+1}-\frac{2}{\Delta x^2}{\eta }_i+\frac{1}{\Delta x^2}{\eta }_{i-1}.$$

$${\left({\partial }_{xx}\eta \right)}^{\text{O4}}={\partial }_{xx}\eta +\mathcal{O}(\Delta x^4),{\left({\partial }_{xx}\eta \right)}^{\text{O2}}={\partial }_{xx}\eta +\mathcal{O}(\Delta x^2)$$

$$\frac{|{\eta }_{i+2}-4{\eta }_{i+1}+6{\eta }_i-4{\eta }_{i-1}+{\eta }_{i-2}|}{12\Delta x^2}=|{\left({\partial }_{xx}\eta \right)}^{\text{O4}}-{\left({\partial }_{xx}\eta \right)}^{\text{O2}}|=\mathcal{O}(\Delta x^2).$$

$${\left({\partial }_{xx}\eta \right)}^{\text{O4}}\approx \frac{15\Delta \eta }{12\Delta x^2},{\left({\partial }_{xx}\eta \right)}^{\text{O2}}\approx \frac{\Delta \eta }{\Delta x^2},$$

$$\frac{|{\eta }_{i+2}-4{\eta }_{i+1}+6{\eta }_i-4{\eta }_{i-1}+{\eta }_{i-2}|}{12\Delta x^2}=|{\left({\partial }_{xx}\eta \right)}^{\text{O4}}-{\left({\partial }_{xx}\eta \right)}^{\text{O2}}|\approx \frac{\Delta \eta }{4\Delta x^2}.$$

$${\sigma }_i=\min (1,\frac{ϵ+\frac{|{\eta }_i-{\eta }_{i-1}|}{\Delta x}+\frac{|{\eta }_i-{\eta }_{i+1}|}{\Delta x}}{ϵ+\frac{|{\eta }_{i+2}-4{\eta }_{i+1}+6{\eta }_i-4{\eta }_{i-1}+{\eta }_{i-2}|}{12\Delta x^2}}),$$

$$\delta_{i}=\frac{1}{4}\hat{\delta}_{i-1}+\frac{1}{2}\hat{\delta}_{i}+\frac{1}{4}\hat{\delta}_{i+1}\,,\quad\hat{\delta}_{i}=\left\{\begin{array}[]{ll}\sigma_{i}&\text{if }\sigma_{i}<0.5\\[10.0pt] 1&\text{otherwise}.\end{array}\right.$$

$$\begin{split}\Delta x&\big{(}\frac{d\mathbf{W}_{i}}{dt}+S_{f_{i}}\big{)}+\dfrac{\mathrm{I}-\mathrm{sign}\big{(}A_{i+1/2}\big{)}}{2}\phi_{i+1/2}+\dfrac{\mathrm{I}+\mathrm{sign}\big{(}A_{i-1/2}\big{)}}{2}\phi_{i-1/2}+\\ \Delta x&\frac{\mu_{i-\frac{1}{2}}}{2}\Big{\{}\frac{1}{3}[\big{(}\frac{d\mathbf{W}}{dt}+S_{f}\big{)}_{i-1}-\big{(}\frac{d\mathbf{W}}{dt}+S_{f}\big{)}_{i}]+[\text{sign}(A)\big{(}\dfrac{d\mathbf{W}}{dt}+S_{f}\big{)}]_{i-\frac{1}{2}}\Big{\}}+\\ \Delta x&\frac{\mu_{i+\frac{1}{2}}}{2}\Big{\{}\frac{1}{3}[\big{(}\frac{d\mathbf{W}}{dt}+S_{f}\big{)}_{i+1}-\big{(}\frac{d\mathbf{W}}{dt}+S_{f}\big{)}_{i}]-[\text{sign}(A)\big{(}\dfrac{d\mathbf{W}}{dt}+S_{f}\big{)}]_{i+\frac{1}{2}}\Big{\}}=\\ f_{b_{i}}&\int\limits_{\Omega}\varphi_{i}D(\mathbf{W}_{\textrm{h}})dx+\sum_{i=0}^{N-1}f_{b_{i+\frac{1}{2}}}\int\limits_{x_{i}}^{x_{i+1}}\left(A(\mathbf{W}_{\textrm{h}})\partial_{x}\varphi_{i}\,\tau_{\textrm{s}}\right)D(\mathbf{W}_{\textrm{h}})dx,\end{split}$$

$$\mu_{i\pm\frac{1}{2}}=\left\{\begin{array}[]{ll}1&\text{if }\;\min(f_{\text{break}_{i}},\,f_{\text{break}_{i\pm 1}})=1\\ \delta_{i\pm\frac{1}{2}}\omega_{i\pm\frac{1}{2}}&\text{otherwise}.\end{array}\right.$$

$${\omega }_i=\text{exp}\left\{-\alpha \Delta x\frac{\underset{j=0,N}{\max }{H}_j-ϵ}{\max (ϵ,\underset{j=i-2,i+2}{\min }{H}_j-ϵ)}\right\},$$

$$\mathcal{E}=\frac{A}{h}.$$

$${\mathcal{E}}_{\text{h}}=\frac{|\eta |}{|h|+ϵ},$$

$$\begin{array}{rl}& |{\partial }_x\eta |\ge \tan \varphi \text{with}\varphi \in [14^{\circ },33^{\circ }]\\ & \frac{{\partial }_t\eta }{\sqrt{g|h|}}>\gamma \text{with}0.35\le \gamma \le \mathrm{0.65.}\end{array}$$

$$\begin{array}{rl}& F{r}_{w_i}^n:=\frac{|{\eta }_i^n-{\eta }_i^{n-1}|}{\Delta t\sqrt{g(|h_i|+ϵ)}}\text{and}\\ & \nabla {\eta }_i^n:=\max \left(\frac{|{\eta }_{i+1}^n-{\eta }_{i-1}^n|}{2\Delta x},\frac{|{\eta }_i^n-{\eta }_{i-1}^n|}{\Delta x},\frac{|{\eta }_{i+1}^n-{\eta }_i^n|}{\Delta x}\right).\end{array}$$

$$F{r}_{w_i}^n>\gamma \text{or}\nabla {\eta }_i^n\ge \tan \varphi ,$$

$$Fr_{b}=\sqrt{\frac{\big{(}2\frac{H_{2}}{H_{1}}+1\big{)}^{2}-1}{8}}.$$

$$F{r}_s=\frac{u_s}{c_b},$$

Full text

We provide a detailed discussion of the implementation of the newly proposed detection method, as well as of two other well known criteria which are used for comparison. An extensive benchmarking on several problems involving different wave phenomena and breaking conditions allows to show the robustness of the numerical method proposed, as well as to assess the advantages and limitations of the different detection methods.

1 Introduction

Wave breaking is a fundamental phenomenon, especially in the nearshore region. It drives wave set-up and set-down, wave run-up, nearshore circulations (e.g. longshore and rip currents), sediment transport, forces on structures, and many other effects. At large scales, an important effect of wave breaking is the energy dissipation, which is the result of many local hydrodynamic phenomena, ultimately leading to dissipation due to small scale viscous effects. In this paper we focus on the wave breaking closure when wave propagation is performed by means of depth averaged Boussinesq models.

Among the main families of closures found in literature we mention the eddy viscosity formulations [52, 38, 25, 42, 14], the breaking roller methods [8, 50], and the hyperbolic or hybrid approaches [6, 47, 44, 23]. Within the first type of approaches, the eddy viscosity one, the macroscopic dissipation of the energy is modelled by means of explicitly introducing a viscous dissipation term in the momentum, or, alternatively, both in the mass and momentum equations. In particular, the eddy viscosity coefficients are defined via physical arguments, or through some auxiliary evolution model (see e.g. [24] and references therein), and, eventually, adjusted by means of numerical experiments. An alternative to the eddy viscosity type modeling are the roller approaches, that account more explicitly for the large scale effects of vorticity on the mean flow. Finally, the hybrid approaches exploit the properties of hyperbolic conservation laws endowed with an entropy, and model large scale effects of wave breaking with the dissipation of the total energy across shocks arising in shallow water simulations or augment/modify Boussinesq-type equations to include, for example, additional terms caused by the presence of the surface rollers (see cf. [41]).

Independently of the chosen closure strategy, some breaking detection criteria is very often required to trigger the onset of wave breaking. In practice it is this criteria that leads to the activation of the closure. As well explained in [33], historically the first detection criteria were based on quantities computed using information over one full phase of the wave. The possibility of performing accurate phase resolved simulations, has led to the necessity of formulating new detection criteria based on local flow features [38, 25, 42, 23]. More specifically, these criteria rely on a wave-by-wave analysis more efficient to program in the context of phase resolved simulations, and providing a physically correct detection of breaking onset and termination, provided that a sufficiently accurate description of the wave profiles is available. In this paper we consider a one-to-one comparison of three different detection criteria, tested on a common solver. The wave propagation is handled by the well known enhanced Boussinesq equations of Madsen and Sørensen [31]. Wave breaking is modelled using a hybrid approach and reverting locally to the non-linear Shallow Water equations. Our focus is devoted to the comparison of numerical implementations of the non-linearity criterion by [45, 46, 47], of the hybrid slope/vertical velocity criterion by [23, 37], and of a convective criterion, which compares the free surface velocity to the celerity of the wave [5]. We propose in particular an implementation in an upwind Finite Element solver, with a novel discontinuity capturing technique, ad-hoc developed for this application. We perform a very thorough benchmarking on flows involving different types of waves and of breaking configurations. The results allow to validate our numerical implementation, and insight on the advantages and limitations of the physical behaviour of the different detection criteria. This paper extends and improves the discussion of [2, 3].

The manuscript is organized as follows. In section §2 we discuss the model equations, while section §3 is devoted to the discussion of the finite element scheme used to solve them. Some details on the detection and capturing of shocks are provided, as they have been developed ad-hoc for this application, and provide an element of novelty. Section §4 encompasses the definition of the breaking detection criteria, and provides a precise description of their numerical implementation. Finally, in section §5 we collect a thorough comparison of the detection criteria on different kinds of waves, along to experimental data. The paper is concluded by a discussion and outlook on future work.

2 Mathematical Model

2.1 Enhanced Boussinesq Equations

The underlying enhanced Boussinesq model used in this work is the well known system by Madsen and Sørensen reading [30]

[TABLE]

Here, following Figure 1, $h$ denotes the depth at still water level, $\eta$ the free surface level, $u$ the depth-averaged velocity, $H=\eta-h\geq 0$ the water depth, $q=Hu$ the volume flux, and $g$ the gravitational acceleration.

The coefficient $C_{f}$ controls the friction, and all dispersive effects are accounted for by the term $\mathcal{D}(\eta,q)$. In particular, the coefficient values are taken as $\beta=1/15$ and $B=\beta+1/3$. This choice allows to match (in absence of friction) a fourth order Padé development of the dispersion relation of the Euler equations [30].

2.2 Wave Breaking Closure: a Hybrid Boussinesq Type Model

To account also for the effects of wave breaking, system (1) must be coupled to some closure model. In the following, we consider a hybrid approach consisting in reverting locally to the non-linear Shallow Water equations, obtained from (1) with $\mathcal{D}=0$. The rationale behind this approach is that in the region of shallow waters, the waves will locally steepen and rapidly turn into shocks, across which the volume flux and the momentum are conserved, while the total energy is dissipated. The interested reader can refer to [6, 45, 23, 44, 24] for further justification.

To implement this closure, we modify the model equations (1) to

[TABLE]

with the definition of $\mathcal{D}(\eta,q)$ as shown in (1). The newly introduced term $f_{\text{break}}$, denotes a so-called flag, that allows to switch from a propagation region ($f_{\text{break}}=1$) to a breaking one ($f_{\text{break}}=0$).

For planar waves (one spatial dimension) there are two key elements to define $f_{\text{break}}$: the detection of breaking onset and termination; the definition of the region in which the breaking closure is activated. This paper focuses mainly on the study of the detection criteria. Note that, although the analysis is performed in the context of hybrid wave breaking closure, detection criteria are an essential element for other approaches as well. The reader might refer to the extensive references provided in [23, 24] for further details. Our results will also be relevant for multi-dimensional flows, which will however require to embed some additional elements related to the direction of wave propagation, and on water level variations in the transversal direction.

3 Numerical Approximation

3.1 Non-linear Upwind Finite Element Discretization

We propose a non-linear variant of the upwind stabilized finite element method discussed in [36]. The method exploits a standard $C^{0}$ finite element approximation on a tessellation of the spatial domain $\Omega=[0,L]$, consisting of $N$ constant intervals, or cells. The cell width is denoted by $\Delta x$. Let us introduce the vector of unknowns $\mathbf{W}=[\eta\,,q]^{T}$, along with

[TABLE]

Denoting by $\varphi_{i}(x)$ (i=0, …,N) the standard linear Lagrange basis functions, the method relies on the expansions

[TABLE]

Introducing the residual

[TABLE]

The discrete equations at a node $i$ are obtained from the following stabilized variational statement

[TABLE]

where $A=\partial F/\partial\mathbf{W}$ denotes the shallow water flux Jacobian, and $\tau_{\textrm{s}}$ a stabilization matrix. For the latter it is well known (cf. [36] and references therein for details), that by defining the matrix as111Matrix absolute values are computed, as usual, by means of an eigen-decomposition. $\tau_{\textrm{s}}=|A|^{-1}\Delta x/2$, the resulting scheme can be recast as the following upwind splitting method with a non-diagonal mass matrix

[TABLE]

where we have introduced the hydrostatic fluctuations $\phi$, with $\phi_{i+1/2}=F_{i+1}-F_{i}+gH_{i+1/2}[0,\,h_{i+1}-h_{i}]$ and $\phi_{i-1/2}$ approximated analogously. The terms $\mathbf{W}_{i\pm 1/2}$ denote cell arithmetic averages.

As for the right hand side terms, neglecting the boundary conditions, which will be discussed later, these are evaluated using the following approximations (cf. definition of $D$ in the beginning of the section and of $\mathcal{D}$ in (1))

[TABLE]

and

[TABLE]

The above expressions are readily evaluated by numerical quadrature in each cell using the linear finite element expansions to interpolate all the quantities with the “h” subscript. Note that, for any quantity $u$ we have $\partial_{x}u_{\textrm{h}}=u_{i+1}-u_{i}/\Delta x$ over a cell $[x_{i},\,x_{i+1}]$. The auxiliary variables $w^{\eta}$ and $w^{q}$ are introduced to handle the higher order derivatives. Their values are obtained by locally inverting the projections

[TABLE]

The interested reader can refer to [36] or more details.

3.2 Shock Capturing, Limiters, and Entropy Fix

In absence of friction and of dispersive effects, if the mass matrices on the left hand side of (5) are lumped, we end up with the first-order upwind flux splitting approximation

[TABLE]

with the short notation $A_{i\pm 1/2}=A(\mathbf{W}_{i\pm 1/2})$. The first-order scheme above is the basis of many discontinuity capturing approximations for shallow water flows (see e.g. [11, 15, 12]). The strategy proposed here is to introduce a cell-wise mass-matrix limiter, that allows to switch from this first order non-oscillatory method across shocks, to the full finite element approximation, which has been shown to be third-order accurate and is appropriate for propagation [36]. Taking into account the explicit form of the mass matrix entries (see [36] for details), the full non-linear approximation can be written as

[TABLE]

For the cell-wise limiter $\delta_{i\pm 1/2}$ we have tested several classical slope/flux limiters. In our case, however, it is very important that numerical dissipation, in correspondence of smooth extrema, is not excessive, as this may hinder the accuracy of the breaking closure. For this reason we have in practice used a limiter based on a smoothness indicator, with the main idea to combine two different approximations of a solution derivative. To fix ideas, let us consider the case of a second-order derivative, but for the first-order derivative one proceeds similarly. We consider a fourth-order approximation

[TABLE]

along with a second-order approximation

[TABLE]

Clearly for a smooth variation of $\eta$ we retrieve

[TABLE]

and hence

[TABLE]

On the other hand, for a discontinuity of finite amplitude $\Delta\eta$ around node $i$ we get

[TABLE]

or, equivalently,

[TABLE]

We thus define the term

[TABLE]

with $\epsilon$ a small parameter used to avoid any division by zero. In the smooth case, while the numerator in the right slot is bounded, the denominator is of $\mathcal{O}(\Delta x^{2})$, and thus the value of $1$ is retained. Across a discontinuity, the numerator gives $\Delta\eta/\Delta x$, while the denominator gives approximately $\Delta\eta/4\Delta x^{2}$, so $\sigma_{i}=4\Delta x\ll 1$. In practice, in every cell we set $\delta_{i\pm 1/2}=\min(\delta_{i},\delta_{i\pm 1})$ where, to avoid abrupt variations of the sensor, we set in every node

[TABLE]

The advantage of this definition is its ability to better preserve smooth extrema, while still capturing strong discontinuities. We have tested this smoothness sensor limiter on several test cases, comparing it to more classical limiters. We report hereafter the behaviours observed in two typical cases. The first consists in the run-up of a solitary wave, which will be discussed in detail in section §5.2. To approximate the solution, we consider the Shallow Water equations with the first-order scheme and test different limiting strategies. In Figure 2, we show the free surface distributions obtained at $t^{\prime}=15$ (cf. §5.2 for more details) for the smoothness sensor, the Superbee and the monotonized centered [28] limiters. Comparing these approximations, one can clearly see that the smoothness limiter improves the preservation of the wave peak.

As a second example, we consider a dam break problem. The results of Figure 3 show that the smoothness sensor allows to capture the strong shock, with a behaviour somewhat close to that provided by the monotonized centered limiter. Note that, to obtain these results, the upwind flux splitting terms $\text{sign}(A)\Delta F$ need to be corrected with an entropy fix, for which we have adapted the techniques detailed in [20, 26, 34] (see [3] for details). These two benchmark problems allow to highlight the advantage of using the proposed smoothness limiter, and all further test cases will be performed relying on this choice. In particular, the complete discretization in its final form reads:

[TABLE]

having denoted by $\mu_{i\pm\frac{1}{2}}$ a cell limiter computed as

[TABLE]

Further, $f_{b_{i}}=f_{\text{break}_{i}}\omega_{i}$, where $\omega(H)$ is the exponential wet/dry cut-off introduced in [35], and computed as

[TABLE]

with $\epsilon$ a mesh-dependent small threshold in practice set to $\epsilon=(\Delta x/L)^{2}$, where $L$ is a characteristic length of the domain, and $\alpha$ is usually set to $\alpha\approx 10$. In practice we set $f_{b_{i\pm 1/2}}=\min(f_{b_{i\pm 1}},f_{b_{i}})$ and similarly for $\omega$.

The resulting scheme can be readily shown to be well-balanced in wet regions. In particular considering the lake at rest state given by $q=0$ and $\eta$ constant across the entire domain, one can show that all the dispersive terms, as well as the hydrostatic fluctuations $\phi_{i\pm 1/2}$ and the friction terms, are identically zero. The scheme thus reduces to $d\mathbf{W}/dt=0$, which shows that this state is indefinitely preserved. More details of the analysis can be found in [36]. In this work we have generalized the validity of this property, to account also for dry areas, adopting the modifications proposed in [11, 15, 35]. Finally, the integration in time of (9) is performed with a non-dissipative implicit Crank-Nicholson method, and wave generation and boundary conditions are handled using the source functions and absorbing layers discussed to some extent in [36]. The interested reader can refer to this reference, and to [3], and references therein, for more details.

Remark 1

Note that the choice of the presented method is not crucial. However, its detailed description allows to underline the strict interplay between the numerics (limiting procedure), and the modelling of wave breaking. While other schemes have been used in combination with a similar wave breaking closure (cf. [23, 19, 24]), such interaction is always necessary. This makes the understanding of the stabilisation and shock-capturing mechanisms critical.

Remark 2

While the numerical approximation has been extensively tested for convergence on the Madsen and Søresen equations in [36], our focus here shall be on the breaking and shock capturing features.

4 Wave Breaking Detection and Implementation

We discuss in this section the detection criteria considered in this paper. The objective is to provide computable definitions for the breaking flag $f_{\text{break}}$ required in the scheme. In particular, we consider three approaches: one based on wave non-linearities, taken from [47, 48]; the second defined on both a slope and a vertical velocity condition, taken from [23]; the third being our implementation of the so-called convective or critical free surface Froude criteria discussed e.g. in [29] (see also [5, 9] and references therein). To our knowledge this is a first attempt at using such a convective criteria in predicting wave breaking in the context of Boussinesq simulations with a breaking closure.

From a physical point of view, it is well known that one of the main parameters ruling wave breaking is the Froude number. However, when looking at wave breaking, several definitions of this quantity are possible, going from the classical local value $Fr=|u|/\sqrt{gH}$, to definitions based on the ratio of the wave amplitude over the depth (cf. [47], or [13] for open channel flows), to other definitions involving wave by wave analysis as the local trough Froude criteria proposed in [33], and many others.

For moving bores, the well known studies by Favre [17] and Treske [49] show that the limiting value of the local Froude number, relative to the bore above which breaking occurs, is about $1.3$. In general, the onset of breaking has been identified to happen in a range of the Froude numbers between $1.3$ and $1.6$ (see [33, 45, 46, 47, 48]). In practice, defining a local Froude number starting from data issued from numerical (or experimental) data is not a trivial task as many parameters are less clearly identified than, e.g., in the experimental set-up of [17] and [49]. In a wave by wave analysis, for example, there is the question of determining the wave’s celerity, the location in which the velocity and Froude condition should be estimated, which velocity should be used, and so on.

The three criteria considered here allow to cover good part of the existing detection methods. To be able to fairly compare them, our objective is to provide approximations to the sought solution, based on the same numerical method, and within the same breaking closure approach. As much as possible of their practical implementation is also discussed, to enable for reproducibility.

4.1 A Criterion Based on Local Non-linearity

The first criterion we have considered is the simplest to set up. It is taken from [45, 46, 47, 48] in which the authors present a detection method based on the non-linearity parameter $\mathcal{E}$, defined as the ratio wave amplitude over still water depth, i.e.

[TABLE]

For the onset of breaking, the value $\mathcal{E}=0.8$ is suggested. This is obtained assuming a transition at a Froude $\approx 1.6$ [33, 45, 48], and considering that the limit value to obtain a stable solitary wave is $\mathcal{E}\approx 0.78$ [22]. Clearly, this criterion is purely local, as it involves no variations (derivatives) of the wave characteristics. Concerning the termination criterion, once $\mathcal{E}>0.8$, the breaking closure is kept active until $\mathcal{E}$ is back within a range between $0.55\div 0.25$, as suggested in [47, 48].

We will refer in the following to this detection method as the Local Criterion. It has the advantage of being simple, but in some cases it lacks generality, as, for example, in presence of wave set-up (or set-down), and in general when the value of $h$ is not anymore relevant to define a reference depth or wave celerity.

In our implementation of this detection method, we have computed in each node the discrete analogue of (10) given by

[TABLE]

with $\epsilon\approx(\Delta x/L)^{2}$, which allows to avoid any division by zero when considering the run-up on sloping beaches. Further, whenever a point in space satisfies $\mathcal{E}>\mathcal{E}_{\text{critical}}$, with $\mathcal{E}_{\text{critical}}=0.8$, we seek the point downstream in the wave direction satisfying $\mathcal{E}=0.3$. The downstream direction is either known a priori, or determined using the technique discussed in section §4.3. In the whole region between the two points we set $f_{\text{break}}=0$.

4.2 Hybrid Slope/Vertical Velocity Criterion

The second criterion we consider was initially proposed in [23]. It combines two criteria: one related to the steepness of the wave [38] and the second involving the vertical displacement of the free surface [7]. In practice, the onset of breaking is determined by the validity of one of the two following conditions

[TABLE]

In [23], the value $\gamma=0.6$ is prescribed as a default. Note that the second condition closely resembles a Froude number condition based on the velocity of the free surface. This condition is inefficient for steady hydraulic jumps, in which case the slope criterion is more effective. For this the choice of the critical angle depends strongly on the breaking type. As in [23] and [44], and unless explicitly mentioned, we have used the default value $\phi=30^{\circ}$. The termination of breaking is based in this case on the Froude condition proposed in [44], which relies on the analysis of moving bores. Following [23], this detection method is referred to as Hybrid Criterion.

The practical implementation is done in four phases as follows.

Pre-flagging. We compute in every node

[TABLE]

If

[TABLE]

we set $f_{\text{break}_{i}}=0$. Adjacent cells involving flagged nodes are clustered, thus obtaining a preliminary definition of the breaking length. 2. 2.

Peak and trough water level. In each cluster of breaking cells we determine the maximum and minimum values of the free surface elevation, $H_{2}$ and $H_{1}$, and the corresponding mesh points $x_{\max}$ and $x_{\min}$, as on Figure 4. 3. 3.

Termination check. We use the criteria proposed in [44] and compute the bore Froude number $Fr_{b}$ defined as

[TABLE]

If $Fr_{b}<Fr_{b_{\text{cr}}}$, we set $f_{\text{break}}=1$ along all the breaking length and move to the next cluster of breaking cells. We set $Fr_{b_{\text{cr}}}=1.3$, as proposed in [23, 44]. 4. 4.

Roller magnitude estimation and final flagging. If $Fr_{b}>Fr_{b_{\text{cr}}}$, we follow again the definitions of [23, 44]: we compute the physical roller length as $l_{r}=2.9(H_{2}-H_{1})$, and we define a numerical roller of length $l_{NSW}=2.5l_{r}$. We set $f_{\text{break}_{i}}=0$ in all points verifying $2x_{i}-(x_{\max}+x_{\min})<\bar{l}_{NSW}$.

4.3 The Convective or Critical Free Surface Froude Criterion

The main idea behind this detection method is that wave breaking occurs when the free surface velocity at the peak of the wave exceeds the celerity of the wave. Therefore, one may consider the Froude number

[TABLE]

with $u_{s}$ the free surface velocity at the wave crest, and $c_{b}$ the celerity. If $Fr>1$ then the wave should break. In the following we will refer to this detection method as to the Physical Criterion.

To implement this method in practice, we must provide computable definitions of both the free surface velocity and the celerity. Concerning the first, we recall the asymptotic development of the horizontal velocity reading [27]

[TABLE]

where $u$ is the depth averaged velocity. The idea is then to evaluate the above expression on the free surface, for $z=\eta=H+h$, and at the wave’s crest, where it is assumed that $\partial_{x}u=0$. Evaluating the above expression, we obtain

[TABLE]

Concerning the celerity of the wave, we combine the two relations

[TABLE]

obtaining

[TABLE]

With these definitions, we proceed to an implementation based on the following four steps.

Pre-flagging. Same as for the hybrid criterion (see section §4.2). 2. 2.

Peak and trough water level. Same as for the hybrid criterion (section §4.2). 3. 3.

Evaluation of $Fr_{s}$. We need to evaluate the free surface velocity in correspondence of $x_{\max}$, as well as the celerity. For the discrete free surface velocity $u_{s_{\textrm{h}}}$, we use (15) in which the second derivative of the velocity is computed first by means of standard second-order centred finite difference formula, and then smoothed, in order to remove some of the numerical noise. The filtering is done using one iteration of a 5-point moving average, using 2 mesh points on either side of $x_{\max}$. We neglect the second term in (15) in all the computations performed, which involve piecewise constant slopes. The implementation of (16), used to evaluate the celerity, is

[TABLE]

with $\epsilon$ added to avoid any division by zero. The critical free surface Froude is then discretely evaluated as

[TABLE] 4. 4.

Roller magnitude estimation and final flagging. If $Fr_{s_{\textrm{h}}}>Fr_{s_{\text{cr}}}$, we proceed as in section §4.2: we define $l_{NSW}=2.5l_{r}$ with $l_{r}=2.9(H_{2}-H_{1})$, and we set $f_{\text{break}_{i}}=0$ in all points verifying $2x_{i}-(x_{\max}+x_{\min})<\bar{l}_{NSW}$. The critical value by default used is $Fr_{s_{\text{cr}}}=1$, but, as we shall see in the numerical experiments, other values have been tested to compensate for the shoaling behaviour of the considered Boussinesq model.

Remark

The determination of the wave celerities necessary for the convective criterion is a key point for future extensions. This issue is common to both numerical and experimental techniques aiming at characterising in detail the wave-field. The filtering approach used here is a first solution, but more involved techniques will be for sure required to implement this approach in a more general context, both concerning the quality of the estimated celerities and the extension to multiple space dimensions.

4.4 Additional Implementation Details

As discussed in detail in [24], the wave breaking closure used here introduces easily oscillations at the interface of the coupling between the Bousisnesq model and the Shallow Water equations. Indeed, in [24], it is shown that this is especially important when a mesh refinement is involved, as, if uncontrolled, it would lead to the blow up of the computation, instead of the sought mesh convergence. In particular, on coarse meshes the amplitude of these perturbations remains very small. However, since there is no clear understanding of the link between the mesh size and numerical spurious effects, we have tried to minimize the sources of instability. To this end, some techniques for the minimization of the oscillations have been already mentioned, as, e.g., the use of filtering to evaluate the second derivatives of the velocity. Furthermore, through the determination of the roller magnitude, as, e.g. within the fourth step of the practical implementation of the physical criterion, we have also avoided the introduction of spurious breaking. Specifically, we have throughout avoided the introduction of rollers with a size smaller than the stencil of our numerical scheme, i.e. for $\bar{l}_{NSW}\leq 4\Delta x$, we set back $f_{\text{break}_{i}}=1$. Similarly, if two breaking regions are too close, they are merged, such that, if the distance between two breakers is less than $4\Delta x$, we set $f_{\text{break}_{i}}=0$ also for the nodes in between.

Finally, we avoid any possible wrong detection, as, e.g., the breaking being detected on the rear side of the wave via the check of the fulfilment of $c_{b_{\textrm{h}}}(x_{\max}-x_{\min})\geq 0$.

5 Numerical Tests and Results

5.1 Wei’s Solitary Wave Shoaling on a Slope

In order to study the physics of the wave breaking, we consider the benchmark problems proposed in Wei et al. ([51]), involving solitary waves shoaling over different constant slopes. These tests are used to investigate the physical behaviour obtained with the different breaking criteria.

5.1.1 Case 1: Slope of 1:35

The set-up is given by a solitary wave of amplitude $A=0.2$m over a spatial domain defined within $[-12,40]$m. The bathymetry, set initially with $h_{0}=1$m depth, has a constant slope 1:35 between $0\leq x\leq 35$m and continues till the end of the domain with a flat slope. The solitary wave starts at the toe of the slope set at $x=0$m. The total considered time is $30$s. The time step is dictated by the CFL number $\nu=0.3$ on a mesh with $\Delta x=0.05$m. Absorbing boundary conditions are set on the left and right ends of the domain using sponge layers $3$m wide (see [36] for details). The wave’s evolution is described as a function of a dimensionless time $t^{\prime}=t\sqrt{g/h_{0}}$.

As in [51], we will focus on four different discrete times: $t_{1}^{\prime}=16.24$, $t_{2}^{\prime}=20.64$, $t_{3}^{\prime}=24.03$ and $t_{4}^{\prime}=25.94$. We visualize the wave profile at these times when using the three different detection criteria: local (Figure 5(a)), hybrid (Figure 5(b)), and physical (Figures 5(c)).

From the figures we can see that all criteria provide the same wave profile at $t_{1}^{\prime}$ and $t_{2}^{\prime}$, at which the wave propagates and starts shoaling, becoming taller and less symmetric. The shoaling continues at $t_{3}^{\prime}$ with the hybrid and physical criteria, but with the local criterion the smooth wave peak has been already flattened, which is the sign of an early onset of wave breaking. This becomes even clearer at time $t_{4}^{\prime}$, at which this criterion provides a wave shape very close to that of a shock. The hybrid criterion fares better, however, the first signs of an early activation of the breaking closure are outlined. Indeed, the wave at $t_{3}^{\prime}$ is still smooth, however we can see the wave height reducing from $t_{4}^{\prime}$ can be seen as the wave profile is already shorter than the previous one, a sign that breaking has taken over the shoaling phase. The physical criterion, instead, shows wave profiles still growing from $t_{3}^{\prime}$ to $t_{4}^{\prime}$, meaning that shoaling is still the dominating phenomenon. As a general remark, compared to the plots in [51], the wave shapes are shorter and less steep. This is related to the fact that the Boussinesq model (1) under-shoals compared both to fully non-linear models and to other weakly non-linear models known to be prone to over-shoaling, as, e.g., the one used in the reference. This is a known phenomenon. The interested reader can refer to [18], where the behaviour of well-known weakly non-linear Boussinesq models, including the one used here, is studied in the transition between weakly and fully non-linear waves. To get further insight into the behaviour of the different solutions, we have singled out the first time of breaking for all the considered detection criteria. In Figure 6, we visualize the wave’s position and shape, as well as the closure region at the first breaking instant computed with the procedure described in section §4. As shown in Figure 6(a), the local criteria provides a very early first onset of wave breaking at $t^{\prime}=23.02$. The figure also shows that the breaking termination condition used within this criteria gives a very short sized roller, which does not reach the foot of the wave. Figure 6(b) provides the same visualization for the hybrid criteria. The onset of breaking is at $t^{\prime}=25.48$, much closer to the reference one. The size of the roller is, in this case, larger, and allows to cover the whole front of the wave, including its foot. This is most likely a better situation for the stability of the criteria, as the steeper parts of the wave are all embedded in the roller. Passing to the physical criterion, we can see in Figure 6(c) that for $Fr_{s_{\text{cr}}}=1$, this method provides a late onset of breaking at $t^{\prime}=28.04$. This is a consequence of the undershoaling properties of the Madsen and Sørensen model, that provide waves with reduced steepness, and thus curvature at the peak. This leads to a late explosion of the free surface velocity $u_{s}$ in (15). The lateness in the onset of the breaking also leads to a considerably taller wave, and, as a consequence, a roller thickness much larger than that observed, e.g., in Figure 6(b). To compensate for the behaviour of the Boussinesq model, several solutions are possible. One could be, for example, to modify (15) in order to compensate for the aforementioned properties of the model, e.g. by replacing the vertical development $u(z)$ by a development on $q(z)$, using ideas discussed in [18]. Here we adopt a simpler solution, and introduce a reduced value of the critical Froude, set to $Fr_{s_{\text{cr}}}=0.75$.

As shown in Figure 6(d) this provides a breaking time $t^{\prime}=25.91$, within $0.1\%$ of the reference, and a slightly reduced roller thickness. The wave profiles obtained in this case are identical to those of Figure 5(c). In [51] a very interesting analysis is reported visualizing the dynamics of the wave in terms of evolution of the peaks celerity versus the particle speed at the peak. This is particularly relevant here, as the physical criterion is actually based on the location where these two curves cross each other. In Figure 7 we compare the temporal evolution of the non-dimensional wave’s celerity $c_{b}^{\prime}=c_{b}/\sqrt{gh_{0}}$ to the non-dimensional particle velocity of the wave’s crest $u_{s}^{\prime}=u_{s}/\sqrt{gh_{0}}$. Note that in the reference the authors report the celerity obtained from $\hat{c}_{w}^{\prime}=\frac{dx_{c}}{dt}/\sqrt{gh_{0}}$, where $x_{c}$ corresponds to the position of the wave’s crest. Here, we use both for $c_{b}$ and for $u_{s}$ the discrete definitions proposed in section §4.3, so that the breaking onset for the physical criterion is exactly represented on the plot.

In figure 7, we also represent with red circles the first activation of breaking (represented as a circle in the reference as well). Even though neither of the underlying models, nor the computation of the celerity and of the particle speed are the same as in [51], the behaviours observed are very close to those of the reference.

As the wave shoals, with all the criteria we obtain an increase in the particle velocity at the crest, and a decreasing celerity at the front. This is expected for shoaling conditions. The results allow to see graphically the early onset of breaking for some of the criteria, which is, of course, natural as the only one actually using the physical quantities involved is the physical criterion. It is interesting to note, however, that after the start of the breaking, all closure models provide a strong increase in the particle velocity at the crest, which is absent in the reference. This is certainly a consequence of the wave profile sharpening, and eventually converging into a shock, and of the use of (15) to compute $u_{s}$. Interestingly, very similar excursion in the free surface velocities in breaking waves can be observed experimentally, see e.g. [32].

5.1.2 Case 2: Slope of 1:15

We consider the same problem set-up, with a solitary wave of amplitude $A=0.3$m, shoaling on a slope 1:15. Following [51], we start looking at the dimensionless time instants $t_{1}^{\prime}=3.23$, $t_{2}^{\prime}=6$, $t_{3}^{\prime}=8.4$ and $t_{4}^{\prime}=11.32$. Differently from the previous case, these instants do not include a breaking point for the reference [51]. In our case, all criteria provide, nevertheless, breaking before $t_{4}^{\prime}$, as shown in Figure 8.

As before, the earliest onset is obtained with the local criterion for $t^{\prime}=7.67$. The hybrid detection method predicts a first breaking at $t^{\prime}=9.32$, while the physical criterion gives breaking onset times of $t^{\prime}=10.99$ and $t^{\prime}=10.12$ with $Fr_{s_{\text{cr}}}=1$ and $Fr_{s_{\text{cr}}}=0.75$ respectively. The latter thus provides the latest breaking onset, and gives wave profiles very close to the ones reported in [51], including the peak within the last wave which has not dissipated yet. These observations are confirmed by the wave profiles at the first breaking instant, reported in Figure 9. As in the previous case, we see that the local criterion gives a smaller roller which does not involve the foot of the wave, while the roller obtained with the hybrid criterion encompasses the whole front of the wave. The physical criterion with $Fr_{s_{\text{cr}}}=1$ provides a taller wave with a roller length which touches the end of the wet area, while with $Fr_{s_{\text{cr}}}=0.75$ a reduced roller is obtained. For both values of $Fr_{s_{\text{cr}}}$ the wave front is sharper than the one obtained with the other two criteria. As before, we compare the time evolution of the crest’s free surface velocity with the wave’s celerity. The results are summarized in Figure 10. The breaking point of the local and hybrid criterion occur very early. The dissipation of the wave is such that the curves never cross. This may be due to the effects of activating the limiter in the shock smoothing of the waves peak. Conversely, the physical criterion provides a similar behaviour as the one observed in the previous case, with a blow-up of the particle velocity at the crest after breaking onset. This is consistent with the sharper wave fronts seen in Figure 9.

5.2 Breaking and Run-up of a Solitary Wave

This is a classical benchmark due to by [43], which adds further complexity to the previous one, and considers the whole transformation of a solitary wave, including the processes of breaking, run-up, and backwash. The set-up consists in a solitary wave of amplitude $A=0.28$m propagating on an initial depth of $h_{0}=1$m and shoaling on a 1:19.85 slope. A sketch of the test is shown in Figure 11.

Computations have been run with a Manning coefficient of $0.01$, a CFL$=0.3$ and $\Delta x=0.1$m, i.e. the same parameters suggested in [23]. The solution is visualized at several discrete instants, corresponding to the beginning of the shoaling, the shoaling, the run-up and the back-wash. Results are displayed w.r.t. the non-dimensional time $t^{\prime}=t\sqrt{g/h_{0}}$, and compared to experimental data. The same problem has been used in section §3 to evaluate the smoothness detection strategy proposed in the paper (see Figure 2).

The results are summarized in Figure 12. Despite the differences observed in the previous case, in this test we can hardly see any distinction between the solutions obtained with the considered detection methods. Overall, each approach shows an excellent agreement with the experimental data. We can see a slight phase advance at $t^{\prime}=15$, which, however, is independent of the detection criterion, as none of the computations has started breaking at this point. This phase difference is assumed to be linked to the celerity of the exact solitary wave for (1) given by (see [36] for details)

[TABLE]

The effect of this phase advance is not visible in the solutions at later times. We underline the nice capturing of the hydraulic jump at $t^{\prime}=45$, which validates the numerical limiting procedure proposed (cf. section §3). All of the breaking criteria manage to capture this feature. The exception to this is the solution obtained with the local criterion, that outlines some weak oscillations, revealing a weak flagging of this region. Last, we would like to remark that for this type of problems a singularity is encountered for the local criterion as $h$ goes to zero when the slope crosses the line corresponding to the initial/still water depth. This singularity is handled by the addition of a small quantity at the denominator, as discussed in sections §4.1 and §4.2. More importantly, the value of $h$ in vicinity and after this point is not relevant any longer for the definition of the physical quantities used in the breaking citeria. In particular, we can note that the local criterion always breaks after this point, as $|\eta|=|H-h|>|h|$ when $h<0$. Similarly, the first condition in the pre-flagging (12) loses its physical meaning, as $\sqrt{gh}$ no longer represents a relevant celerity. For cases as the one in this section, this does not seem to affect the results. This is due to the fact that the three criteria end up detecting breaking in the run-up phase, and close to the backwash region where the depth is very small. These features are very well reproduced by the Shallow Water equations. However, as we will see later, for periodic waves, and especially in presence of strong set-up, this may be a source of inaccuracies.

5.3 Periodic Wave over a Submerged Bar

In this benchmark we consider monochromatic waves propagating over a submerged bar. According to [40], this test originates from Dingemans [16] who verified Delft Hydraulics numerical model HISWA, and was then repeated by Beji and Battjes later [4].

The test consists in a periodic wave with amplitude $A=0.027$m and period $T=2.525$s. The wave is generated at $X_{0}=22$m over an initial water depth of $h_{0}=0.4$m. The bathymetry has a fore slope of 1:20, and a rear one of 1:10. The length of the plateau is of $2$m, and the corresponding depth of $0.1$m. Experimental data in four gauges (cf. Figure 13) are available for comparison. The problem is solved on the domain $[0,60]$m, with $\Delta x=0.04$m, and absorbing conditions on both ends. Periodic waves are generated using the source method of [36]. This is a difficult case, and following [23], we modify the values of section §4 for the detection criteria. We set $\mathcal{E}_{\text{critical}}=0.6$ for the local criterion, and $\gamma=0.3$, $\phi=30\degree$ for the hybrid and physical ones. For the latter, we will use only $Fr_{s_{\text{cr}}}=0.75$.

For $Fr_{s_{\text{cr}}}=1$ no breaking is observed at all. Moreover, reflections due to the boundary conditions may induce a phase shift in the results, affecting the comparison with the experimental data [21]. For this reason, as in [23] we have focused our attention on the first generated waves, before any interaction with the boundary conditions on the right end of the domain has taken place. In practice, the first gauge is used to re-phase the numerical signals, and the resulting phase calibration is applied to all gauges. The resulting water elevation time series in the three remaining gauges are reported in Figure 14.

The comparison shows that the hybrid criterion is the only one capturing the correct magnitude and phase of the waves. The local criterion essentially fails to capture the onset of breaking, even with the reduced value of $\mathcal{E}_{\text{critical}}=0.6$. This leads to wrong amplitudes as well as to a shift in the phase of parts of the signal. The physical criterion gives acceptable results, managing to capture the phase of the experimental data. The amplitudes remain slightly higher than those of the data and obtained with the hybrid detection. This may be the consequence of a later onset of breaking, as well as of some intermittency in the breaking process, reducing the overall dissipation.

To better understand the differences in the time series, we have plotted the wave profiles at time instants corresponding to the first and last breaking instants on the plateau. Note that these differ for the three criteria. We also have included the wave profile for all the solutions at the same intermediate time. The resulting free surface distributions are displayed in Figure 15.

If we set $t=0$ in correspondence of the first breaking of the hybrid criterion, the physical criterion breaks roughly at $0.18\,T$, while the local detection gives an onset at $t=0.345\,T$, with $T$ the incoming wave period. The delay is clearly visible in the left pictures in Figure 15. The plots also show that the local criterion still gives a very small roller, while the hybrid and physical criterion gives a breaking region of comparable size. The central pictures, referring to $t=0.365\,T$, show that the local criterion has already deactivated the flagging, while both the hybrid and physical criteria are still breaking. Before the wave leaves the plateau, the local criterion activates again at $t\approx 0.82\,T$, probably due to the decrease in bathymetry. This sudden flash in breaking has no effect on the wave shape. The hybrid and the physical criteria keep the flagging activated until $t\approx 0.77\,T$s and $t\approx 0.70\,T$, respectively. One must however also consider that in general, at least in our implementation, the lack of memory in the flagging introduces a degree of intermittency. This is observed more for the physical criterion than for the hybrid one. The cause may be the noise in the evaluation of the derivatives defining $u_{s}$ (cf, section §4.3), and the result is a weaker dissipation leading to the slightly higher amplitudes of Figure 14.

Mesh convergence study

We use this test to evaluate the mesh dependence of the closure and detection criteria. Note that a very thorough investigation of this issue has already been presented in [24]. The reference, as well as our study shows that the type of breaking closure based on the coupling of Boussinesq-type equations with the Non-linear Shallow Water model and the approximation of the rollers by bores, introduces instabilities, that lead to a blow up, as the mesh gets refined. This behavior corresponds most likely to a mismatch between the energy transported in the dispersive field which is not entirely transmitted at the interface. Smoothing the coupling numerically only delays this blow up. Note that similar results, albeit with different schemes, are reported in [39, 19, 24].

The grid convergence study reported in this section allows to include in the grid convergence analysis the effect of using different detection criteria. For the scheme and model used here, grid convergence for smooth, non-breaking, cases can be already found in [36, 19]. We thus repeat the computations of the submerged bar by halving the mesh size from $\Delta x=0.04$m, used in the previous section, to $\Delta x=0.02$m and $\Delta x=0.01$m. Only, the local criterion has allowed to compute the whole solution with an additional refinement of $\Delta x=0.005$m.

On Figure 16 we report the results for the water levels in gauge $G_{3}$, but the trends observed in the other gauges are similar. In the figure, the first row reports the results obtained with the local criterion, the second those of the physical detection criterion, and the last the results obtained with the hybrid criterion. For the last two criteria our results are very similar to those of [24], and show the appearance on the finer meshes of unphysical high frequency oscillations superimposed to the trend of the coarse mesh solution. This phenomenon is a not as visible for the results obtained with the local criterion, which only flashes at the end of the upward slope. Compared to the weakly nonlinear model analyzed in [24] the oscillations obtained here seem slightly weaker in both in frequency and amplitude. We believe this to be a consequence of the shoaling properties of the Madsen-Sørensen equations used here which have a tendency to under-shoal, compared to the Nwogu-model of the reference, and to give lower and longer waves, especially in the nonlinear range [18]. Concerning the detection criteria, the main conclusions we can draw are that the onset and growth of the instabilities seem to be more visible with the hybrid detection which for this case is the one breaking more strongly. The phenomenon is still very intense with the physical criterion, which still detects wave breaking on most of the plateau, and much less present with the local detection criterion which breaks very weakly. Indeed, the local criterion, for the mesh $\Delta x=0.005$m appears to have the same shape as a Boussinesq-type solution, i.e. without any breaking applied at all, as can be seen from Figure 17.

Finally, as a last comparison, in Figure 18 we show how the physical and hybrid criterion start failing with the very first breaking instances, when the mesh gets too fine, as with here $\Delta x=0.005$m.

5.4 Breaking and set-up of Monochromatic Waves on a Shore

As a last benchmark, we consider the shoaling and breaking of monochromatic waves over a shore, based on the experiments of Hansen and Svendsen [10]. This is an interesting case, as the set-up is induced by the breaking near the shore. The waves considered are generated on a still water depth $h_{0}=0.36$m and then propagate on a slope of 1:32.26. As in [23] the grid size used is $\Delta x=0.025$m. Wave generation, as well as the sponge layers on the left end of the domain, are set up following [36].

Case 051041. For this case waves have a period $T=2$s and amplitude $A=0.018$m, corresponding to an Ursell number $U=4.8077$. We visualize in Figure 19 the wave profiles at the first and last breaking instants (for a fixed wave) detected by each criteria.

As before, the breaking region is highlighted by vertical lines. As in the other cases, the onset of breaking is predicted slightly later by the physical criterion when $Fr_{s_{\text{cr}}}=1$. Differently than in the other benchmarks, here for $Fr_{s_{\text{cr}}}=0.75$ the onset is predicted earlier than with all other methods. The termination is also shifted accordingly for the two cases. We can also see again that the roller size computed by the local criterion is much shorter. For the latter, we can see a single vertical line at $x\approx 11.5$, corresponding to the zero of $h$, where breaking is turned on until the end of the domain. As discussed at the end of section §5.3 this is a consequence of the use of $h$ in the definition of a reference depth in the breaking detection. In Figure 20, we compare the spatial distribution of the wave height and the mean water level (set-up) with the experimental data of [10]. The results are qualitatively close to those reported in [23]. The prediction of the wave heights is clearly affected from the prediction of the breaking onset. Even though under-predicting the breaking wave height, the physical criterion with $Fr_{s_{\text{cr}}}=1$ gives the best result for this case in terms of position of the breaking onset and slope of the height after the breaking point. All the other criteria, including the physical one for $Fr_{s_{\text{cr}}}=0.75$, provide a very early onset of breaking, with an under-prediction of the breaking wave height, and a milder slope after. Some of the height underprediction may be also related to the under-shoaling characteristics of the Madsen and Sørensen model (1). Looking at the wave set-up, reported on the right in the same figure, we can clearly see the effects of the use of $h$ in the local criterion. An unrealistically high set-up is observed. As shown in [6] this can be related to the amount of dissipation injected in the region where $h<0$ which is entirely flagged as breaking. This has no effects on the results obtained with the other criteria, providing a correct set-up.

Case 031041. In this case waves have a period $T=3.333333$s and amplitude $A=0.0215$m, corresponding to an Ursell number $U=17.5588$.

The results are summarized in Figures 21, and 22, reporting respectively the wave profiles at breaking onset/termination for a given wave, and the mean water levels and set-up. We can observe a behaviour similar to the one discussed in the previous tests. The local criterion gives very thin rollers, and flags as breaking occurs within the entire region where the bathymetry is above the zero of $h$. The onset is predicted slightly later for the physical criterion if $Fr_{s_{\text{cr}}}=1$, but the difference is less important in this case. Concerning the mean water level, the best is again the one obtained with the physical criterion and $Fr_{s_{\text{cr}}}=1$, the other methods, including the physical criterion for a lower value of the critical Froude, are affected by the early prediction of breaking onset.

The same drawbacks of the local criterion can be observed in the wave set-up, which is too large for this detection method, at least with the implementation discussed in section §4. A good prediction of this quantity is instead provided by the hybrid and physical criteria.

6 Conclusions and Outlook

We have discussed the numerical implementation and the extensive benchmarking of three breaking detection criteria in a finite element Boussinesq model. Breaking closure is handled by means of the hybrid method, that reverts the enhanced Boussinesq equations to the hydrostatic Shallow Water system, and models wave breaking fronts as shocks. The detection criteria tested involve an original implementation in a numerical scheme of a more physical detection based on the comparison of the free surface particle velocity at the peak of the wave, with the wave celerity. In the case of the convective criterion, the lowering of the $Fr_{s_{\text{cr}}}$ to $0.75$ has been set empirically, to account for the under-shoaling of the Boussinesq model [18]. The thorough benchmarking performed shows an advantage in using detection methods based on more sound physical arguments. In its current implementation the physical criterion is comparable to the hybrid slope/vertical velocity method proposed in [23], which is already quite satisfactory in a wide range of cases. It also allows to obtain improved results in some cases, as the Hansen and Svendsen experiments of section §5.4. The results also show some limitations related to the adequacy of the detection criteria to the propagation model. More specifically, the classical asymptotic formula used for the vertical velocity profile may be less well adapted for the type of Boussinesq model used than for, e.g., such as the standard Peregrine equations of the Beji and Nadaoka [18]. This is clearly an aspect which could be refined and investigated in the future, also including fully non-linear models. Further, the switch between the detection criteria undergoes some smoothing but this will not guarantee the stability of the coupling as $\Delta x\rightarrow 0$. This has been already investigated in [19, 24] for the Green-Naghdi and Nwogu Boussinesq equations. The grid convergence study performed here confirms this observation also for the Madsen and Sorensen Model. The study also allows to give some highlight on the impact of the nature of flagging on this aspect. In particular more than the intermittency of the flagging, the results suggest that the stronger the duration of the breaking process the earlier the breakdown. When convergence is obtained, as in the case of the local criterion in our test, the converged solution is essentially the non-breaking Boussinesq one, which we have shown quantitatively. The breakdown of the solution as the mesh is refined is a flaw for this closure, which deserves a specific investigation.

Another element which deserves further study is the definition of the front celerity, which becomes less trivial in multiple space dimensions. This may require to combine the ideas explored here with techniques used in image processing, as well as, in wind waves analysis. The main challenge will be to retain the locality of the detection, thus avoiding non-local operations as, e.g., Fourier transforms. Other elements more specifically linked to the numerical method also deserve some further refinement. It is the case for example for the approximation of the second derivatives necessary in the vertical velocity development. These may require special attention in the multidimensional case, especially on unstructured grids. The use of methods based on high-order polynomial expansions (see [1]) may help, but some form of filtering will certainly still be necessary.