<figure><img src="image_1.png" /></figure>

HAL
open science

# Exploring a conservative staggered scheme for Boussinesq-type equations: Insights into numerical diffusion, dispersion, and wave-breaking

Fatima-Zahra Mihami, Volker Roeber

► To cite this version:

Fatima-Zahra Mihami, Volker Roeber. Exploring a conservative staggered scheme for Boussinesq-type equations: Insights into numerical diffusion, dispersion, and wave-breaking. Coastal Engineering, 2025, pp.104880. 10.1016/j.coastaleng.2025.104880. hal-05296160

HAL Id: hal-05296160

https://hal.science/hal-05296160v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Exploring a conservative staggered scheme for Boussinesq-type equations:
Insights into numerical diffusion, dispersion, and wave-breaking

Fatima-Zahra Mihami b,a, ID,* , Volker Roeber a, ID

## ARTICLE INFO

Dataset link: https://github.com/fzmihami/WavePropError

**Keywords:**

Boussinesq-type modeling

Numerical dispersion and diffusion

Wave breaking

Turbulent kinetic energy (TKE)

Finite difference and finite volume methods

Grid sensitivity

Conservative staggered scheme

Numerical stability

## ABSTRACT

Accurate and efficient modeling of coastal wave transformation, particularly under wave-breaking conditions, remains a major challenge for Boussinesq-type models. To address this, we introduce and validate a conservative staggered-grid scheme to discretize a set of weakly nonlinear Boussinesq-type equations. The presented approach revisits the staggered finite-difference strategy by ensuring a momentum-conserving solution designed to enhance numerical stability and improve shock-capturing properties. The scheme's performance is assessed through a series of numerical tests involving monochromatic and spectral linear wave propagation. These tests demonstrate that the conservative staggered scheme is much less sensitive to grid resolution, resulting in approximately one order of magnitude lower numerical diffusion in contrast to the well-established HLLC scheme, while maintaining comparable dispersive accuracy despite using a lower-order spatial reconstruction. Additionally, the scheme introduces a slight negative phase error that compensates for the positive dispersion error inherent in the underlying equations, resulting in improved overall phase accuracy relative to the HLLC scheme. Beyond linear wave propagation, the numerical approach is validated against standard benchmark tests with solitary and spectral breaking waves. In these highly non-linear cases, coupling the conservative staggered scheme with a turbulent kinetic energy (TKE)-based eddy viscosity model yields localized and physically consistent dissipation while preserving the dispersive characteristics of the solution. Compared to conventional hybrid breaking approaches, the TKE-based closure provides enhanced stability, reduced grid sensitivity, and a more accurate representation of energy dissipation during wave breaking. These results underscore the potential of the conservative staggered scheme as an efficient and robust framework for computing complex coastal and nearshore wave processes.

## 1. Introduction

Ongoing sea-level rise driven by climate change has significantly increased the vulnerability of nearshore regions to wave-induced flooding and related coastal hazards. Accurately resolving wave dynamics in these areas is therefore critical, particularly where interactions with complex coastal topography occur over short spatial scales. In this context, the Boussinesq-type approximation, derived from Joseph Boussinesq's 19th-century theory of weakly dispersive waves (Boussinesq, 1872), has become a widely adopted and efficient computational framework for nearshore wave modeling. It captures essential wave transformation processes with high accuracy while offering a substantial reduction in computational cost compared to full three-dimensional flow solvers. Various formulations of Boussinesq-type equations have been developed to improve their accuracy and extend their applicability to a broader range of nearshore hydrodynamic conditions. While

the classical formulation derived by Peregrine (1967) retains only the lowest-order nonlinear and dispersive terms, more advanced formulations have included higher-order components, as shown in the works of Dingemans (1973), Madsen and Schäffer (1998), Karambas and Memos (2009). Alternatively, other approaches have focused on optimizing the dispersion and shoaling characteristics by introducing free parameters into the equations, without resorting to higher-order derivatives, as proposed by Madsen and Sørensen (1992) and Nwogu (1993). Such equations, in both weakly and fully nonlinear forms, now serve as the mathematical foundation for many operational nearshore wave models due to their robustness and efficiency (e.g., FUNWAVE-TVD (Shi et al., 2011), BOSZ (Roeber and Cheung, 2012), COULWAVE (Kim et al., 2009), TUCWave (Kazolea and Delis, 2013), MIKE 21 (Sørensen et al., 2004)).

<footnote>ª Université de Pau et des Pays de l'Adour, E2S-UPPA, SIAME, Anglet, France</footnote>

<footnote>ᵇ University of Hawai'i at Mānoa, Department of Oceanography, Honolulu, HI, USA</footnote>

<footnote>* Corresponding author at: University of Hawai'i at Mānoa, Department of Oceanography, Honolulu, HI, USA.
E-mail addresses: fzmihami@hawaii.edu (F.-Z. Mihami), volker.roeber@univ-pau.fr (V. Roeber).</footnote>
The numerical implementation of Boussinesq-type models has been subject to many trends. Earlier solutions were primarily based on Finite Difference (FD) methods (Abbott et al., 1978; Madsen and Sørensen, 1992; Wei and Kirby, 1995; Beji and Nadaoka, 1996). This approach initially gained traction due to its simplicity and straightforward application. Two main FD strategies emerged, including implicit staggered schemes based on Alternating Direction Implicit (ADI) methods (e.g., Abbott et al. (1978), Madsen and Sørensen (1992)), and collocated grid approaches combined with fourth-order Adams-Bashforth-Moulton (ABM) time integration (e.g., Wei and Kirby (1995), Lynett et al. (2002)). However, FD methods, particularly those using unstaggered grids, often produced numerical noise and spurious oscillations near shorelines and in rapidly varying regions such as breaking wave crests (Kirby, 2016). To address these issues, many models incorporated ad hoc treatments, including filtering and artificial smoothing techniques, to stabilize the solutions (e.g., Madsen et al. (1997), Lynett et al. (2002)).

Recently, a significant number of Boussinesq-type models have adopted a hybrid approach that combines Finite Volume (FV) and FD methods (Erduran et al., 2005; Tonelli and Petti, 2009; Roeber et al., 2010; Kazolea and Delis, 2013; Shi et al., 2012). The decision to transition to FV calculations was driven by the robustness and shock-capturing capabilities of these conservative methods, especially in handling sudden changes in the flow condition, such as the case for wave-breaking and wet-dry transitions. The implementation of the hybrid FV-FD method involves rewriting the dispersive governing equation in a conservative form where the hyperbolic component of the equations is solved using Godunov-type methods paired with Approximate Riemann solvers (e.g., HLL (Harten et al., 1983), HLLC (Toro and Chakraborty, 1994), Roe (Roe, 1997), Osher (Engquist and Osher, 1981)). These solvers are designed to preserve the hyperbolicity of the equations to allow for the formation of discontinuities in the numerical solution. In contrast, the dispersive terms are computed using centered FD approximations on a collocated grid. To enhance spatial accuracy, various high-order reconstruction techniques along with Total Variation Diminishing (TVD) slope limiters have been proposed (e.g., Jiang and Shu (1996), Kim and Kim (2005), Choi et al. (2018)). In recent years, several Boussinesq-type models have used high-order flux reconstruction methods to reduce the dissipative effects of their numerical solution. The fourth-order MUSCL-TVD scheme, introduced by Yamamoto and Daiguji (1993), has found widespread application in numerous Boussinesq implementations, serving as a favored method for approximating the dispersive conservative equations (e.g., Erduran et al. (2005), Tonelli and Petti (2009), Shi et al. (2012)). Another successful approach, relying on the Weighted Essentially Non-Oscillatory (WENO) scheme described by Jiang and Shu (1996), has been used to solve the weakly nonlinear equations of Nwogu (1993), as shown in Zhou et al. (2016). Regarding the time integration, the choice of the time-stepping algorithm for the Boussinesq-type models has been evenly split between the Strong Stability-Preserving (SSP) Runge-Kutta schemes (Roeber and Cheung, 2012; Shi et al., 2012; Dutykh et al., 2011; Cienfuegos et al., 2006) and ABM schemes (Erduran et al., 2005; Tonelli and Petti, 2009).

Another numerical approach that has received little attention within the Boussinesq community, despite its widespread use in non-hydrostatic Shallow Water model implementation (e.g., SWASH (Zijlema et al., 2011), NEOWAVE (Yamazaki et al., 2012) and Xbeach (Roelvink et al., 2018)), is the conservative staggered scheme. This method, originating from the work of Stelling and Duinmeijer (2003), Stelling and Zijlema (2003), combines the efficiency of the staggered FD approximation with the shock-capturing properties achieved through a momentum-conserving discretization. Such numerical techniques have been successfully used to solve problems involving high Froude number flows, such as bores, hydraulic jumps, and the inundation of dry land (e.g., Zhou and Stansby (1999), Stelling and Duinmeijer (2003), Madsen et al. (2005), Doyen and Gunawan (2014), Yamazaki et al. (2009)). The

scheme employs momentum-conservative approximations that satisfy the Rankine-Hugoniot jump condition at a discrete level, allowing for valid solutions in rapidly varying flows and effectively computing the dynamics of traveling bores associated with wave breaking in the surf zone (see Zijlema (2019)). Moreover, the use of second-order upwinding in the momentum equations effectively eliminates short-wave instabilities, ensuring model stability without the need for artificial filtering (Zijlema et al., 2011). For these reasons, the conservative staggered scheme appears as an attractive alternative to the traditional Godunov-type methods that rely on computationally expensive approximate Riemann solvers and high-order flux reconstruction.

Over the past 15 years, the use of shock-capturing schemes in Boussinesq-type models has gained significant importance, especially for handling wave-breaking. These recent developments have been closely tied to the implementation of a technique that is conventionally referred to as the “hybrid method”, in which the dispersive terms are deactivated when a certain criterion for wave breaking is met. This reduces the dispersive model momentarily to a cell-wise hyperbolic Shallow Water model, allowing for the breaking wavefront to evolve as a moving hydraulic jump. First introduced by Tonelli and Petti (2009), this wave-breaking method relies on the conservation of mass and momentum across discontinuities and thus goes hand in hand with the use of the FV methods that solve the hyperbolic portion of the Boussinesq-type equations. Although the implementation of such a method is rather straightforward, the stability and accuracy of the results greatly depend on the breaking criteria, which determine the transition from a dispersive to a hydrostatic state. Many breaking onset criteria have been proposed based on theoretical studies, numerical simulations, laboratory experiments, or field observations (e.g., Barthelemy et al. (2018), Tissier et al. (2012), Roeber et al. (2010), Kazolea and Ricchiuto (2018)). The hybrid wave-breaking method has also proven to be very sensitive to the grid resolution, as noted by Kazolea and Ricchiuto (2018). The abrupt transition between dispersive and non-dispersive solutions can introduce spurious oscillations and contaminate the solution, potentially leading to instabilities as the grid is refined (Filippini et al., 2016). The stability of the solution herein often relies to some degree on numerical diffusion to overcome the intrinsically incompatible coupling of the dispersive/non-dispersive states. In contrast, other wave-breaking closures, particularly those employing the eddy viscosity approach, have shown less sensitivity to grid size. This approach introduces additional dissipation to the momentum equations to mimic the energy dissipation that arises from turbulence and mixing (Zelt, 1991). While the eddy viscosity source term is relatively straightforward, assessing its magnitude remains challenging. The term should accurately reproduce the correct amount of energy dissipation and must act locally at the breaking wavefront. Kennedy et al. (2000) introduced a formulation for the eddy viscosity term, which is expressed as a function of the time derivative of the free surface elevation. A more sophisticated wave-breaking closure stems from the Prandtl-Kolmogorov one-equation approach (i.e., Kolmogorov (1941), Prandtl et al. (1945)) that computes the temporal and spatial variation of turbulent kinetic energy (TKE) in an additional transport equation and subsequently makes use of the Boussinesq concept by relating TKE to eddy viscosity. Although the TKE-based eddy viscosity closure has a larger computational cost compared to the hybrid technique, it has shown lower sensitivity to mesh size and greater stability and robustness for both weakly and fully nonlinear Boussinesq-type equations (see Kazolea and Ricchiuto (2018)).

Within the framework of Boussinesq-type modeling, the choice of the numerical method for the discretization can have a significant impact on the validity and computational efficiency of the numerical solution. Aiming for shock-capturing capabilities, several recent Boussinesq-type models have been relying on Godunov's type methods for the solution of the hyperbolic part of the equations. Though these schemes provide enough numerical dissipation to mask high-frequency spurious oscillations, the stability is gained at the cost of
wave-damping (Kirby, 2020). Even with high-order reconstruction, Godunov's type methods inherently introduce diffusion and consequently require a high-resolution grid to maintain wave amplitude during propagation. A grid sensitivity study with the Boussinesq model FUNWAVE-TVD shows the effects of spatial resolution on energy loss during wave propagation (Torres et al., 2022). At least 60-grid points per wavelength are recommended for the setup of the model. Another limitation of this numerical method is the high computational cost required by the approximate Riemann solvers in combination with the high-order flux reconstruction that must be computed at each grid interface across the domain. Yuan et al., 2020 showed that in their FUNWAVE-TVD implementation, the computation of the flux with the Riemann solver takes around 42.9% of the total model computation time on the CPU. Hence, using a simplified flux approximation may significantly reduce the overall computational cost of a model.

In this paper, we propose an alternative numerical framework for the solution of the Boussinesq-type equations. The approach is based on the conservative staggered scheme proposed by Stelling and Duinmeijer (2003) and enhanced to achieve second-order accuracy in both space and time. The weakly dispersive Nwogu equations (i.e., Nwogu (1993)) are chosen as the mathematical basis for the new implementation, alongside a TKE-based eddy viscosity method to deal with wave breaking. Here, we contrast this new numerical approach to the standard hybrid FV-FD method commonly used in many Boussinesq-type models (e.g., Shi et al. (2011), Roeber and Cheung (2012), Kim et al. (2009), Kazolea and Delis (2013)). A series of one-dimensional linear wave tests is conducted to evaluate numerical dissipation and dispersion errors, followed by more realistic one- and two-dimensional benchmark cases to assess the accuracy of the new implementation in capturing nonlinear processes, including wave breaking. The remainder of this paper is organized as follows. Section 2 presents the governing equations of Nwogu (1993) in conserved variable form and details the discretization for the staggered grid. Section 3 provides a numerical analysis of dispersion and diffusion errors, using linear wave solutions for simplicity and comparing the scheme against the well-established HLLC method. Section 4 describes the validation of the conservative staggered scheme through one- and two-dimensional benchmark cases involving wave-breaking and nonlinear wave processes. Finally, Section 5 presents the discussion and conclusion.

## 2. Numerical methods

### 2.1. Governing equations

The dispersive Boussinesq-type equations, derived by Nwogu (1993), consist of continuity and momentum equations. These equations describe the motion of the free surface and its two orthogonal velocity components in the 2D horizontal plane while maintaining a pseudo-3D structure by incorporating additional terms to account for non-hydrostatic pressure effects. The system can be written in conservative vector form:

$$U_t + F(U)_x + G(U)_y + P(U) + S(U) = 0 \quad (1)$$

$U$ is the conservative variable vector, $F(U)$ and $G(U)$ are the advective acceleration vectors, $P(U)$ is the depth-integrated hydrostatic pressure and $S(U)$ is the source term which includes dispersion terms:

$$U = \begin{bmatrix} H \\ P \\ Q \end{bmatrix} \qquad F(U) = \begin{bmatrix} H u \\ H u^2 \\ H u v \end{bmatrix} \qquad G(U) = \begin{bmatrix} H v \\ H u v \\ H v^2 \end{bmatrix} \qquad P(U) = \begin{bmatrix} 0 \\ g H \eta_x \\ g H \eta_y \end{bmatrix} \quad (2)$$

$$S(U) = \begin{bmatrix} \psi_c \\ u \psi_c - H_t \psi_u^n + H[\psi_u^c]_t \\ v \psi_c - H_t \psi_v^n + H[\psi_v^c]_t \end{bmatrix}$$

where:

$$P = H \left( u + \frac{z_{\alpha}^2}{2} u_{xx} + z_{\alpha} (hu)_{xx} \right) \quad (3)$$

$$Q = H \left( v + \frac{z_{\alpha}^2}{2} v_{yy} + z_{\alpha} (hv)_{yy} \right) \quad (4)$$

$h$ denotes the still water depth, assumed to be stationary (i.e., $h_t = 0$). $\eta$ is the elevation of the free surface and $H = \eta + h$ is the total depth of the water. $(u, v)$ is the horizontal flow velocity at $z_{\alpha}$. The reference depth $z_{\alpha}$ is the vertical position within the still water column at which the velocity is evaluated. Nwogu (1993) defines this level as a fixed proportion of the local water depth with a position near mid-depth: $z_{\alpha} = \sigma h = -0.531h$.

The dispersive contribution to the continuity equations is represented by the term $\psi_c$, defined as:

$$\psi_c = \left[ \left( \frac{z_{\alpha}^2}{2} - \frac{h^2}{6} \right) h (u_{xx} + v_{xy}) + \left( z_{\alpha} + \frac{h}{2} \right) h ((hu)_{xx} + (hv)_{xy}) \right]_x + \left[ \left( \frac{z_{\alpha}^2}{2} - \frac{h^2}{6} \right) h (u_{xy} + v_{yy}) + \left( z_{\alpha} + \frac{h}{2} \right) h ((hu)_{xy} + (hv)_{yy}) \right]_y \quad (5)$$

The momentum dispersion variables are divided into cross- and normal-derivative terms:

$$\psi_u = \psi_u^n + \psi_u^c \quad \text{where} \quad \begin{cases} \psi_u^n = \frac{z_{\alpha}^2}{2} u_{xx} + z_{\alpha} (hu)_{xx} \\ \psi_u^c = \frac{z_{\alpha}^2}{2} v_{xy} + z_{\alpha} (hv)_{xy} \end{cases} \quad (6)$$

$$\psi_v = \psi_v^n + \psi_v^c \quad \text{where} \quad \begin{cases} \psi_v^n = \frac{z_{\alpha}^2}{2} v_{yy} + z_{\alpha} (hv)_{yy} \\ \psi_v^c = \frac{z_{\alpha}^2}{2} u_{xy} + z_{\alpha} (hu)_{xy} \end{cases} \quad (7)$$

A key point of the conservative formulation used in this paper is that it avoids splitting the free surface gradient into an artificial flux gradient and a source term, which account for bed slope. As a result, the surface gradient is calculated independently from the numerical flux, eliminating the need for additional treatment to ensure the well-balance of the numerical solution (see Mihami et al. (2022)).

### 2.2. Spatial discretization

Eq. (1) forms a hyperbolic equation with the dispersion serving as a source term. In this section, we describe the spatial discretization of the equations that relies on a conservative staggered scheme. The state variables are approximated on a staggered C-grid, where the free surface elevation $\eta$ is defined at the cell center, and the depth-integrated velocity $(u, v)$ is estimated at the cell interface. Conservative variables $[H, P, Q]^T$ are defined in a similar way (see Fig. 1).

To describe the numerical discretization, we consider a 2D rectangular computational domain with a uniform grid spacing of $\Delta x$ and $\Delta y$ in the $x$ and $y$ directions, respectively. The variables stored at the cell center are expressed as $x_{i,j}$, where $i$ and $j$ are the spatial indices in the $x$ and $y$ directions. The variables stored at the cell interface are indicated by $x_{i\pm\frac{1}{2},j}$ or $x_{i,j\pm\frac{1}{2}}$ in the $x$- and $y$-directions, respectively. Time stepping is based on discrete, non-uniform time levels $t^n$, where the adaptive time step is defined as

$$\Delta t^n = t^{n+1} - t^n.$$

Each time level is computed recursively as

$$t^n = \sum_{k=0}^{n-1} \Delta t^k, \quad \text{with} \quad t^0 = 0.$$

The value of a variable $a$ at time level $t^n$ is denoted by $a^n$. The water depth $H$ is evaluated at each time level $t^n$, while the depth-averaged velocities $u$ and $v$ are evaluated at intermediate times $t^{n+\frac{1}{2}} = \frac{1}{2}(t^n+t^{n+1})$.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Definition sketch: Nwogu's equations state variables.</figcaption></figure>

This temporal staggering facilitates consistent second-order accuracy in both space and time.

The spatial discretization of the governing equation (i.e., Eq. (1)) can be outlined as follows:

$$ \left[ \frac{\partial U}{\partial t} \right]_{i,j} + \frac{1}{\Delta x} \left( F_{i+\frac{1}{2},j} - F_{i-\frac{1}{2},j} \right) + \frac{1}{\Delta y} \left( G_{i,j+\frac{1}{2}} - G_{i,j-\frac{1}{2}} \right) + P_{i,j} + S_{i,j} = 0, \quad (8) $$

Here, we define the spatial discretization of the state variable $\mathcal{V}$ and conservative variables $U$ as:

$$ U_{i,j} = \begin{bmatrix} H_{i,j} \\ P_{i+\frac{1}{2},j} \\ Q_{i,j+\frac{1}{2}} \end{bmatrix}, \quad \mathcal{V}_{i,j} = \begin{bmatrix} H_{i,j} \\ u_{i+\frac{1}{2},j} \\ v_{i,j+\frac{1}{2}} \end{bmatrix} \quad (9) $$

The other terms in the governing equation follow the same index notation.

### Convection terms

In Eq. (8), the flux contributions $F_{i+\frac{1}{2},j}$ and $G_{i,j+\frac{1}{2}}$, represent the discrete convection terms. These are computed with a second-order upwind scheme. In this case, the choice of upwind variables and the upwinding criterion are crucial to satisfying the Rankine-Hugoniot jump condition at a discrete level (Zijlema, 2019). As long as the discrete fluxes satisfy the Rankine-Hugoniot relation, discontinuities such as bores and shocks are preserved without spurious oscillations, which makes the scheme inherently shock-capturing (see Stelling and Duinmeijer (2003)). We reformulate convection terms as a multiplication of an unwinded variable $\mathcal{V}$ and an upwinding vector.

$$ F_{i+\frac{1}{2},j} = \mathcal{V}_{i+\frac{1}{2},j} \cdot \begin{bmatrix} u \\ Hu \\ Hu \end{bmatrix}_{i+\frac{1}{2},j} \qquad \mathcal{G}_{i,j+\frac{1}{2}} = \mathcal{V}_{i,j+\frac{1}{2}} \cdot \begin{bmatrix} v \\ H\nu \\ H\upsilon \end{bmatrix}_{i,j+\frac{1}{2}} \quad (10) $$

The state variable $\mathcal{V}$ is approximated at the cell interface using a second-order upwind scheme:

$$ \mathcal{V}_{i+\frac{1}{2},j} = \begin{cases} \mathcal{V}_{i,j} + \frac{1}{2}\psi \left( \mathbf{r}_{i+\frac{1}{2},j}^{+} \right) \left( \mathcal{V}_{i,j} - \mathcal{V}_{i-1,j} \right), & \text{if } [u, Hu, Hu]_{i+\frac{1}{2},j}^{T} \ge 0 \\ \mathcal{V}_{i+1,j} + \frac{1}{2}\psi \left( \mathbf{r}_{i+\frac{1}{2},j}^{-} \right) \left( \mathcal{V}_{i+1,j} - \mathcal{V}_{i+2,j} \right), & \text{if } [u, Hu, Hu]_{i+\frac{1}{2},j}^{T} < 0 \end{cases} \quad (11) $$

$$ \mathcal{V}_{i,j+\frac{1}{2}} = \begin{cases} \mathcal{V}_{i,j} + \frac{1}{2}\psi \left( \mathbf{r}_{i,j+\frac{1}{2}}^{+} \right) \left( \mathcal{V}_{i,j} - \mathcal{V}_{i,j-1} \right), & \text{if } [v, H\upsilon, H\nu]_{i,j+\frac{1}{2}}^{T} \ge 0 \\ \mathcal{V}_{i,j+1} + \frac{1}{2}\psi \left( \mathbf{r}_{i,j+\frac{1}{2}}^{-} \right) \left( \mathcal{V}_{i,j+1} - \mathcal{V}_{i,j+2} \right), & \text{if } [v, H\upsilon, H\nu]_{i,j+\frac{1}{2}}^{T} < 0 \end{cases} \quad (12) $$

where $\mathbf{r}$ is the gradient of the flow variable $\mathcal{V}$ approximated at the cell interface:

$$ \mathbf{r}_{i+\frac{1}{2},j}^{+} = \frac{\mathcal{V}_{i+1,j} - \mathcal{V}_{i,j}}{\mathcal{V}_{i,j} - \mathcal{V}_{i-1,j}}, \qquad \mathbf{r}_{i+\frac{1}{2},j}^{-} = \frac{\mathcal{V}_{i+1,j} - \mathcal{V}_{i,j}}{\mathcal{V}_{i+2,j} - \mathcal{V}_{i+1,j}} \quad (13) $$

$$ \mathbf{r}_{i,j+\frac{1}{2}}^{+} = \frac{\mathcal{V}_{i,j+1} - \mathcal{V}_{i,j}}{\mathcal{V}_{i,j} - \mathcal{V}_{i,j-1}}, \qquad \mathbf{r}_{i,j+\frac{1}{2}}^{-} = \frac{\mathcal{V}_{i,j+1} - \mathcal{V}_{i,j}}{\mathcal{V}_{i,j+2} - \mathcal{V}_{i,j+1}} \quad (14) $$

$\psi(\mathbf{r})$ is the slope limiter function, which locally reduces the solution from second to first order at locations with sharp gradients. In this work, a MinMod slope limiter is used:

$$ \phi(\mathbf{r}) = \max(0, \min(1, \mathbf{r})) \quad (15) $$

### Hydrostatic pressure term

In Eq. (8), the pressure-gradient source term is denoted by $P_{i,j}$, which corresponds to the hydrostatic pressure contribution. In the context of depth-integrated equations, such as the Boussinesq-type and Shallow Water equations, the hydrostatic pressure is expressed as the gradient of the free surface: $P_{i,j} = [0, gH\eta_x, gH\eta_y]^T$. In the traditional Godunov-type methods, this term is typically divided into a flux, computed using approximate Riemann solvers, and a source term that accounts for bathymetry variations. This is done to maintain the hyperbolicity of the equations. Here we avoid the splitting, and we compute the term with a second-order central FD:

$$ [gH\eta_x]_{i+\frac{1}{2},j} = g \left( \frac{H_{i+1,j} + H_{i,j}}{2} \right) \left( \frac{\eta_{i+1,j} - \eta_{i,j}}{\Delta x} \right) \quad (16) $$

$$ [gH\eta_y]_{i,j+\frac{1}{2}} = g \left( \frac{H_{i,j+1} + H_{i,j}}{2} \right) \left( \frac{\eta_{i,j+1} - \eta_{i,j}}{\Delta y} \right) \quad (17) $$

Together with the convection terms, the hydrostatic pressure gradient contributes to the shock-capturing property of the staggered finite-difference scheme, since both are discretized in a manner that satisfies the Rankine-Hugoniot condition (see Zijlema (2019)).

### Dispersion terms

In Eq. (8), the dispersive contribution is denoted by $S_{i,j}$. These terms are discretized using FD approximations, similar to the well-established
hybrid FV-FD approach (e.g., Roeber and Cheung (2012), Shi et al. (2011)). However, thanks to the spatial staggering of the variables this discretization results in a more compact FD stencil and better accuracy compared to the collocated grid, as shown in Fornberg (1998).

The dispersion term in the continuity equation can be approximated as follows:

$$[\psi_c]_{ij} = \left[ \left( \frac{z_a^2}{2} - \frac{h^2}{6} \right) h(u_{xx} + v_{xy}) + \left(z_a + \frac{h}{2}\right) h((hu)_{xx} + (hv)_{xy}) \right]_{i+\frac{1}{2},j} - \left[ \left( \frac{z_a^2}{2} - \frac{h^2}{6} \right) h(u_{xx} + v_{xy}) + \left(z_a + \frac{h}{2}\right) h((hu)_{xx} + (hv)_{xy}) \right]_{i-\frac{1}{2},j} + \left[ \left( \frac{z_a^2}{2} - \frac{h^2}{6} \right) h(u_{xy} + v_{yy}) + \left(z_a + \frac{h}{2}\right) h((hu)_{xy} + (hv)_{yy}) \right]_{i,j+\frac{1}{2}} - \left[ \left( \frac{z_a^2}{2} - \frac{h^2}{6} \right) h(u_{xy} + v_{yy}) + \left(z_a + \frac{h}{2}\right) h((hu)_{xy} + (hv)_{yy}) \right]_{i,j-\frac{1}{2}} \quad (18)$$

The second-order space derivatives in Eq. (18) are discretized using second-order FD approximations. Similarly, the dispersion terms in the momentum equations are computed at the cell interface with:

$$[\psi_a]_{i+\frac{1}{2},j} = \frac{z_a^2}{2} \left( [u_{xx}]_{i+\frac{1}{2},j} + [v_{xy}]_{i+\frac{1}{2},j} \right) + z_a \left( [(hu)_{xx}]_{i+\frac{1}{2},j} + [(hv)_{xy}]_{i+\frac{1}{2},j} \right)$$

$$[\psi_v]_{i,j+\frac{1}{2}} = \frac{z_a^2}{2} \left( [u_{xy}]_{i,j+\frac{1}{2}} + [v_{yy}]_{i,j+\frac{1}{2}} \right) + z_a \left( [(hu)_{xy}]_{i,j+\frac{1}{2}} + [(hv)_{yy}]_{i,j+\frac{1}{2}} \right)$$

It should be noted that the source terms for each momentum equation must be estimated at the corresponding cell interface. All source terms in the x-momentum equation are computed at $x_{i+\frac{1}{2},j}$, while the source terms in the y-momentum equation are calculated at $x_{i,j+\frac{1}{2}}$.

Since the continuity equation is solved first, the temporal derivative $H_t$ in the momentum equations can be computed explicitly using the previous and updated values of $H$:

$$[H_t]_{i+\frac{1}{2},j}^{n+\frac{1}{2}} = \frac{[H_{i+\frac{1}{2},j}]^{n+1} - [H_{i+\frac{1}{2},j}]^n}{\Delta t^n} \qquad [H_t]_{i,j+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{[H_{i,j+\frac{1}{2}}]^{n+1} - [H_{i,j+\frac{1}{2}}]^n}{\Delta t^n} \quad (19)$$

Similar to the approach adopted in Shi et al. (2012) and Roeber and Cheung (2012), in order to simplify the system of equations that arise from the time-derivative dispersion terms, all cross-derivatives are moved to the right-hand side of the governing equation and are treated as source terms. The cross-terms $\psi_u^c$ and $\psi_v^c$ are saved at each time step, and the time derivative is then computed with a backward difference approximation:

$$[\psi_u^c]_t]_{i+\frac{1}{2},j}^{n+\frac{1}{2}} = \frac{[\psi_u^c]_{i+\frac{1}{2},j}^{n+\frac{1}{2}} - [\psi_u^c]_{i+\frac{1}{2},j}^{n-\frac{1}{2}}}{\Delta t^{n-1}} \qquad [\psi_v^c]_t]_{i,j+\frac{1}{2}}^{n+\frac{1}{2}} = \frac{[\psi_v^c]_{i,j+\frac{1}{2}}^{n+\frac{1}{2}} - [\psi_v^c]_{i,j+\frac{1}{2}}^{n-\frac{1}{2}}}{\Delta t^{n-1}} \quad (20)$$

At the initial time step, the backward difference approximation requires velocity information from a previous step, which is not available. To resolve this, we assume that the time derivative of the dispersive cross terms (i.e., $\psi_u^c$ and $\psi_v^c$) is initially zero. This is physically justified, as all test cases begin from a state of rest with zero initial velocity. As a result, dispersive terms, dependent on velocity gradients, are inactive at $t = 0$ and only become active as the velocity field evolves over time.

## 2.3. Time-integration

To match the second-order accuracy in space, a combination of the Leapfrog scheme and the Runge-Kutta (RK) predictor-corrector method

(i.e., Gottlieb et al. (2001)) is used. The variable $U$ is calculated using a two-stage time integration with an intermediate solution obtained by the prediction step. In the first step, the equations are solved with only the advection terms and the dispersion terms on the right-hand side:

$$U^* = U^n - \Delta t F(U^n) - \Delta t G(U^n) - \Delta t S(U^n) \quad (21)$$

where:

$$U^n = \begin{bmatrix} H^n \\ P^{n+\frac{1}{2}} \\ Q^{n+\frac{1}{2}} \end{bmatrix} \quad (22)$$

This leads to a predictor solution of first-order accuracy for the convective and dispersion terms. In the second step, the surface gradient terms are added to the momentum equations, and the predicted variables are corrected to full second-order accuracy in time by:

$$U^{n+1} = \frac{1}{2}(U^n + U^*) - \frac{\Delta t}{2}[F(U^*) - G(U^*) - S(U^*)] - \Delta t P(U^{n+1}) \quad (23)$$

It is important to note that the predictor step of the Runge-Kutta time integration should only include the convective acceleration and dispersion terms. The free surface gradient $P$ attains second-order accuracy through Leapfrog time integration combined with the staggering of the flow variables in time. In addition, the use of the TVD Runge-Kutta method for the approximation of the convective acceleration has been shown to enhance the stability of the numerical model. In contrast, non-hydrostatic shallow water models such as SWASH and Xbeach have used a second-order MacCormack method to achieve second-order accuracy for the advection terms (see Zijlema et al. (2011), Roelvink et al. (2018)), while other nonhydrostatic models maintain only first-order accuracy for these terms (e.g., Yamazaki et al., 2011).

The explicit time integration of the hyperbolic equations (Eq. (1)) computes the conservative variables. Since the source terms depend on the primitive variables $u$ and $v$, these variables must be computed at each time step. The variables $P$ and $Q$ can be approximated using a second-order accurate central difference scheme:

$$[P]_{i+\frac{1}{2},j} = H_{i+\frac{1}{2},j} \left( u_{i+\frac{1}{2},j} + \frac{z_a^2}{2} \left[ \frac{u_{i+\frac{1}{2},j} - 2u_{i+\frac{1}{2},j} + u_{i-\frac{1}{2},j}}{\Delta x^2} \right] + z_a \left[ \frac{(hu)_{i+\frac{1}{2},j} - 2(hu)_{i+\frac{1}{2},j} + (hu)_{i-\frac{1}{2},j}}{\Delta x^2} \right] \right)$$

$$[Q]_{i,j+\frac{1}{2}} = H_{i,j+\frac{1}{2}} \left( v_{i,j+\frac{1}{2}} + \frac{z_a^2}{2} \left[ \frac{v_{i,j+\frac{1}{2},j} - 2v_{i,j+\frac{1}{2}} + v_{i,j-\frac{1}{2}}}{\Delta y^2} \right] + z_a \left[ \frac{(hv)_{i,j+\frac{1}{2}} - 2(hv)_{i,j+\frac{1}{2}} + (hv)_{i,j-\frac{1}{2}}}{\Delta y^2} \right] \right)$$

The cross-derivative terms are approximated as source terms in the momentum equations ($\psi_u^c$ and $\psi_v^c$), resulting in the decoupling of the x- and y-systems and instead forming independent tri-diagonal systems. The systems are then solved using the Thomas algorithm at each time step twice (predictor and corrector steps) to calculate the velocity variables $u$ and $v$.

The time stepping is constrained by the Courant-Friedrichs-Lewy (CFL) condition and depends on the grid size and the maximum wave speed:

$$Cr = \Delta t \max \left( \frac{|u|_{i+\frac{1}{2},j}^{n+\frac{1}{2}}| + \sqrt{g H_{i+\frac{1}{2},j}^n}}{\Delta x}, \frac{|v|_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}| + \sqrt{g H_{i,j+\frac{1}{2}}^n}}{\Delta y} \right) < 1 \quad (24)$$

All following numerical tests use a Courant number of $Cr = 0.5$.

The numerical scheme includes an embedded wetting-drying algorithm that requires no ad-hoc treatment. It is well-balanced and preserves the positivity of the water depth under the CFL condition,
allowing run-up and inundation to emerge naturally from the solution (see Gunawan (2015)). To avoid excessively high velocity values near extremely shallow depths, momentum updates (i.e., P and Q) are suppressed when the local depth falls below a threshold $H_{min}$, and the velocities are set to zero (see Toro (2001)). All test cases presented here are based on $H_{min} = 10^{-4}$ m.

## 2.4. Wave breaking

Wave breaking, a key mechanism of energy dissipation in the surf zone, can be modeled using a turbulence kinetic energy (TKE)-based eddy viscosity closure. This method has been shown to be both computationally efficient and to exhibit good convergence properties compared to alternative formulations (Kazolea and Ricchiuto, 2018). In this framework, wave breaking is represented by a time- and space-dependent source term introduced into the momentum equations:

$$R_b = \begin{bmatrix} \left[ \frac{v_t H u_x}{x} \right]_x + \frac{1}{2} \left[ v_t H (u_y + v_x) \right]_y \\ \left[ \frac{v_t H v_y}{y} \right]_y + \frac{1}{2} \left[ v_t H (u_y + v_x) \right]_x \end{bmatrix} = \begin{bmatrix} 0 \\ R_x \\ R_y \end{bmatrix} \quad (25)$$

The discretization of the terms $R_x$ and $R_y$, along with the grid indexing strategy, is provided in Appendix A.1. The eddy viscosity $v_t$ plays a central role in modeling wave breaking-induced turbulence, as it governs the momentum dissipation in the surf zone. It is modeled as a function of the turbulent kinetic energy $k$, based on classical turbulence theory. Kolmogorov (1941) and independently Prandtl et al. (1945) proposed an approximation that relates the eddy viscosity to $k$ and a mixing length scale $l_t$:

$$v_t = C_v \sqrt{k} \cdot \ell_t. \quad (26)$$

where, $C_v$ is a constant coefficient that is set to 0.55 according to many authors (Zhang and Wu, 2011; Kazolea and Ricchiuto, 2018; Pope, 2000). The mixing length $l_t$ is set to the local water depth $h$, a common assumption in depth-integrated wave models (Kazolea and Ricchiuto, 2018). Although this represents a simplification, it has proven to yield robust and consistent results across a wide range of wave conditions without requiring case-specific calibration. Accordingly, we adopt $\ell_t = h$ uniformly in all simulations presented in this work. Pope (2000) proposed a transport equation for $k$ to determine its spatial and temporal variation. The eddy viscosity term (Eq. (26)) is then computed and applied to the entire computational domain. Notably, the eddy viscosity can have nonzero values even in regions where the wave-breaking criterion is not locally satisfied. The term can be seen as referring to the whitewater region, which is intense and very turbulent at the location of wave breaking, but quickly settles down. The governing equation for $k$ is given by:

$$k_t = -\mathcal{A} - \mathcal{E} + \mathcal{P} + \mathcal{D} \quad (27)$$

where $\mathcal{A}$, $\mathcal{E}$, $\mathcal{P}$, $\mathcal{D}$ denote advection, elimination/destruction, production, and diffusion of the turbulent kinetic energy. These terms have been described in detail in Nwogu (1996), Kazolea and Ricchiuto (2018) and involve mean flow quantities and several empirical coefficients.

The advection term $\mathcal{A}$ is defined as:

$$\mathcal{A} = uk_x + vk_y \quad (28)$$

The TKE elimination/destruction terms are denoted as:

$$\mathcal{E} = C_d \frac{k^{3/2}}{\ell_t} \quad (29)$$

where, following Nwogu (1996), Zhang and Wu (2011), Kazolea and Ricchiuto (2018), we set this coefficient to: $C_d = C_v^3$. In this study, the coefficients $C_d$ and $C_v$ are fixed for all numerical tests. However, their values can be optimized for specific cases or wave conditions to improve model performance.

The TKE diffusion term is small and given by

$$\mathcal{D} = \nabla \cdot (v_t \nabla k) \quad (30)$$

The final term, $\mathcal{D}$, represents the production of turbulent kinetic energy (TKE). Following the approach of Nwogu (1996), this term is expressed as:

$$\mathcal{D} = B \frac{\ell_t^2}{\sqrt{C_d}} [(u_z|_\eta)^2 + (v_z|_\eta)^2]^{3/2} \quad (31)$$

where $u_z|_\eta$ and $v_z|_\eta$ are the vertical gradients of the horizontal velocities at the free surface. These variables can be approximated using the irrotationality condition ($w_x = u_z$, $w_y = v_z$), which yields:

$$\begin{aligned} u_z|_\eta &= -\eta [u_{xx} + v_{xy}] - [(hu)_{xx} + (hv)_{xy}] \\ v_z|_\eta &= -\eta [u_{xy} + v_{yy}] - [(hu)_{xy} + (hv)_{yy}] \end{aligned} \quad (32)$$

As TKE production occurs exclusively in cells with wave breaking, it is multiplied by a binary flag $B$, which equals 1 when the breaking criterion is satisfied and 0 otherwise. In this work, wave breaking is triggered when the Froude number at the free-surface exceeds unity:

$$Fr|_\eta = \frac{\sqrt{(u|_\eta)^2 + (v|_\eta)^2}}{\sqrt{gh}} > 1 \quad (33)$$

This criterion marks the onset of supercritical flow, i.e., when particle velocity exceeds the wave celerity. While the shallow-water estimate $C = \sqrt{gh}$ is an approximation under dispersive conditions, it remains a widely used and reasonable choice (Varing et al., 2021).

It should be noted that the overall solution is not very sensitive to the choice of the threshold. For example, Bondehagen et al. (2024) used $Fr|_\eta > 0.85$ as originally analyzed by Barthelemy et al. (2018) as the point, beyond which wave breaking inevitably occurs over the course of a wave period.

Once the wave breaking criterion is flagged, the vertical gradients of the horizontal velocities in Eq. (31) are the primary contributors to the magnitude of $\mathcal{D}$. Energetic breaking, such as that seen in plunging breakers, increases the vertical gradient of horizontal velocities under the wave crest, resulting in an increase in $u_z$ and $v_z$.

## 3. Sensitivity of the staggered scheme to dispersive and diffusive error

In the context of Boussinesq-type wave modeling, the choice of numerical scheme can significantly affect the stability, accuracy, and computational efficiency of the numerical solution, especially when dealing with wave propagation, breaking, and run-up. The discrete solution inherent in numerical schemes can introduce both diffusive and dispersive errors, which can compound the errors already present in the Boussinesq-type equations. While this error can be reduced by refining the grid resolution or using high-order numerical methods, doing so results in a significant increase in computational cost, which can be detrimental for large-scale simulations. This challenge is especially pronounced in long-distance propagation scenarios involving periodic waves, such as swells. The continuous sequence of crests and troughs can amplify the numerical error over time, leading to cumulative effects that distort the wave characteristics.

In this section, we examine the impact of numerical errors on the solution of periodic wave propagation over a flat bottom. We contrast the proposed scheme with a Riemann-based Godunov-type scheme (Harten, Lax, van Leer Contact (HLLC), Fraccarollo and Toro (1995)) commonly used in Boussinesq-type modeling. In this paper, we focus on the numerical solution of Nwogu's equations, however, the findings can be generalized to other Boussinesq-type formulations. The main difference between the Godunov-type schemes and the presented scheme lies in the solution of the hyperbolic Shallow Water equations. The dispersion terms, on the other hand, are approximated in the same way for both the Godunov-type and the presented scheme, using Central
Finite Differences. A key distinction between the two approaches is how they handle the hydrostatic pressure term gHηx. As previously discussed, the present scheme avoids splitting this term into a flux gradient and a source term, whereas the Godunov-type schemes require this splitting to preserve the hyperbolicity of the equation to achieve shock-capturing calculations.

The source code developed for this analysis is available on GitHub as WavePropError (Mihami, 2025), along with detailed documentation and illustrative examples. The repository also includes comprehensive descriptions of the Harten–Lax–van Leer Contact (HLLC) scheme, as well as other Riemann-based and Riemann-free schemes with various flux reconstruction methods.

### 3.1. Monochromatic wave propagation

To analyze the dispersive and diffusive errors, we consider the propagation of a monochromatic wave, characterized by a single frequency and amplitude, over a flat bottom. A monochromatic wave is particularly useful for this analysis because it provides a clear and simple scenario where the influence of numerical errors can be distinctly observed. By examining how the numerical scheme handles the wave's amplitude and phase as it propagates, we can assess the dispersion error, which affects the wave celerity and phase shift, as well as the diffusive error, which impacts the magnitude of the wave damping.

The numerical results are compared to an exact reference solution derived from the linearized Nwogu equations. Specifically, we use the analytical form:

$$ \eta(x, t) = A \cos(kx - \omega t), \quad \text{with } \omega = kC_{\text{Nwogu}} \quad (34) $$

where k is the wave number, ω is the angular frequency, and C_Nwogu is the linear phase speed obtained from the dispersion relation:

$$ C_{\text{Nwogu}} = \left\{ g h \cdot \frac{1 - \left( \alpha + \frac{1}{3} \right) (k h)^2}{1 - \alpha (k h)^2} \right\}^{1/2} \quad (35) $$

with $\alpha = \frac{\sigma^2}{2} + \sigma = -0.39$ (see Nwogu (1993)).

In the following test, we compare the present second-order numerical scheme to the well-established HLLC scheme, which employs fifth-order MLP reconstruction (Kim and Kim, 2005) and third-order Runge–Kutta time integration. It is worth noting that other Godunov-type schemes, such as HLL (Harten et al., 1983) and Central Upwind (Kurganov and Petrova, 2007), have been implemented and all yield identical results to the HLLC scheme for linear wave propagation over flat topography. Proving that for a smooth linear wave and constant water depth, these different numerical approaches collapse into the same solution. Moreover, the numerical diffusion appears to be more sensitive to the order of flux reconstruction rather than the Godunov-type scheme (additional examples can be found in the repository documentation (Mihami, 2025)).

It should be noted that prior to this analysis the implementation of both numerical schemes were verified through a solitary-wave propagation test (Appendix A.2). Unlike periodic monochromatic waves, the solitary wave provides a smooth long-wave solution with a single peak, making it less prone to cumulative dispersive or diffusive errors. This limiting case offers a robust baseline check of the solvers, whereas the sinusoidal wave tests presented here examine how numerical errors accumulate over multiple cycles from recurring extrema.

In a long flat channel, we generate and propagate a T₀ = 4 s period monochromatic wave. The water depth is d = 1 m and the wave amplitude is a₀ = 5 × 10⁻³ m, placing the wave in the linear regime based on the well-known Le Méhauté (1976))'s graph and it is revised version by Zhao et al. (2024). The waves are generated inside the computational domain by introducing a source term to the continuity equation, following the approach proposed by Wei et al. (1999). This approach is analogous to adding and subtracting mass. To ensure a

smooth transition of the wave input, the source term is applied over a range of grid points following a Gaussian distribution. The waves are then absorbed at the boundaries using sponge-layer implementation, as described in Larsen and Dancy (1983), with the sponge layer width at both boundaries set to the wavelength λ₀. We run the test with the two numerical schemes with varying grid resolutions (λ₀/Δx = 15, 30, 45). The free surface snapshot at t = 40T₀ is shown in Fig. 2.

As expected, the numerical error increases with a coarser grid (15 grid points per wavelength) for both methods. While both schemes exhibit dispersive error that alters the wavelength during propagation, the HLLC scheme displays a more pronounced diffusive error compared to the conservative staggered scheme, resulting in a significant wave damping after just a few cycles. With a finer grid (i.e., 45 points per wavelength), the difference between the two solutions becomes minimal, and both converge to the exact solution. It is important to note that numerical errors are cumulative, meaning that the error becomes more significant with longer propagation. However, the tendencies of both dispersive and diffusive errors remain consistent for both schemes. To better investigate the error tendencies, we plotted the time series of the free surface at x = 25λ₀ computed with both schemes and compared the results to the exact solution (Fig. 3.) The time series plot provides a more detailed view of the dispersive error. It is important to note that the exact solution shown here represents the linear Nwogu solution, not the more complete Airy wave theory solution. In practical applications, any dispersive error from Nwogu's equation will either cancel or compound the existing numerical error.

The conservative staggered scheme introduces a negative dispersive error, leading to an underestimation of the wave celerity. This is evident in the time-series plots, where the peaks and troughs of the wave occur later than they should. In contrast, the HLLC scheme has a positive dispersive error, which causes an overestimation of the celerity and the wavelength. While, increasing the grid resolution reduces the dispersive error and phase-shift, the overall trend of the error remains consistent across different grid sizes and wave characteristics. This behavior is further confirmed by additional results for a shorter wave case (kd = π) presented in Appendix A.3, which show the same tendencies: negative phase dispersion error with the staggered scheme, positive with HLLC, and negligible diffusive error for the staggered scheme. The only differences are in the magnitude of the phase shift and amplitude decay. This trend can also be compared to the trend of the dispersive error in Nwogu's equations, where the phase speed error is generally positive for shorter waves (see Nwogu (1993)). This implies that the negative dispersive error from the conservative scheme may counterbalance the positive error in Nwogu's equations, while the positive dispersive error from the HLLC scheme could amplify it.

To quantify the numerical error, we record the computed free surface time series at each grid point between x ∈ [15λ₀, 20λ₀]. The wave amplitude aᵢ at each grid point i is extracted using the Hilbert transform: the analytic signal is constructed from the discrete numerical time series, and its instantaneous envelope is averaged to estimate the local amplitude. The reference amplitude a₀ is taken as the target input wave amplitude. The local phase shift φᵢ between the numerical and exact solutions is obtained by computing the time lag corresponding to the maximum of the cross-correlation between the two time-series. The source code developed for this analysis is openly available at WavePropError (Mihami, 2025). The Root Mean Square Error (RMSE) for both amplitude and phase is then computed using the following formulas:

$$ \text{RMSE}_a = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (a_i - a_0)^2} \qquad \text{RMSE}_\phi = \sqrt{\frac{1}{N} \sum_{i=1}^{N} \phi_i^2} \quad (36) $$

where N is the number of grid points between 15λ₀ and 25λ₀, a₀ is the exact wave amplitude, aᵢ is the computed amplitude, and φᵢ is the phase shift of the numerical solution when compared to the exact linear Nwogu's solution.
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Free surface Snapshot for monochromatic wave propagation with $kd = 0.5$ and $a_0/d = 5 \times 10^{-3}$ at $t = 40T_0$: A comparison between the present scheme (second-order in space and time) and the HLLC scheme (third-order time integration and fifth-order reconstruction) for different grid sizes.</figcaption></figure>

To facilitate comparison across different wave scales and conditions, we introduce the scaled RMSE, defined as the relative amplitude error $RMSE_a/a_0$, and the normalized phase error $RMSE_\phi/T_0$, where $T_0$ is the target wave period. These non-dimensionalized metrics allow us to assess the error independently of the absolute wave size or time scale.

The numerical error for the linear wave propagation is influenced by the grid size and wave characteristics, namely, frequency dispersion parameter (i.e., water depth divided by the wave-length) and the non-linearity parameter (i.e., wave height divided by wavelength), when these parameters are fixed, the scaled numerical errors: $RMSE_a/a_0$ and $RMSE_\phi/T_0$, are expected to stay constant across different problem scales. To verify this, we consider a linear monochromatic wave propagation with the following parameters from Le Méhauté (1976):

$\epsilon_0 = \frac{2a_0}{gT_0^2} = 4 \times 10^{-5}$ and $\mu_0 = \frac{d}{gT_0^2} = 0.02$, which corresponds to $\lambda_0 = 6d$. We run four tests, varying the water depth $d$ to 1 m, 10 m, 100 m, and 1000 m, while keeping the values of $\epsilon_0$ and $\mu_0$ constant. For each run we choose a grid size such as $\lambda_0/\Delta x = 25$ and we use the following numerical schemes:

1. Conservative Staggered scheme with second-order accuracy in both time and space.

2. HLLC scheme with 5th-order flux reconstruction and 3rd-order time integration.

3. HLLC scheme with 3rd-order flux reconstruction and 3rd-order time integration.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Free surface time-series for monochromatic wave propagation with $kd = 0.5$ and $a_0/d = 5 \times 10^{-3}$ at $x = 25\lambda_0$. A comparison between the present scheme (second-order in space and time) and the HLLC scheme (third-order time integration and fifth-order reconstruction) for different grid sizes.</figcaption></figure>

4. HLLC scheme with 2nd-order flux reconstruction and 3rd-order time integration.

The HLLC scheme has been combined with 3rd-order RK time integration to ensure stable results for higher-order flux reconstruction, mainly 5th-order. The computed dispersive and diffusive errors for different scales is shown in Fig. 4. The plot shows that for different water depth both the diffusive and dispersive errors remain approximately constant, provided the wave proprieties, $\epsilon$ and $\mu$ are kept constant. We can notice some slight deviation in the phase error for most of the numerical methods, which can be attributed to small errors from the wave generation and absorption implementations. In the case of HLLC scheme with second-order reconstruction, the scaled RMSE diffusive error, $RMSE_a/a_0$, is around 0.84, indicating a loss of 84% of the wave amplitude after only 20 cycles of propagation. The dissipation error for the HLLC scheme is reduced when using higher-order reconstruction,

dropping to 34% and 19% of the original wave amplitude for third-order and fifth-order, respectively. For the same grid size ($\lambda_0/\Delta x = 25$), the dissipation error with the conservative staggered scheme is approximately 0.2%, which is two orders of magnitude smaller than that of the HLLC scheme, which uses a fifth-order flux reconstruction and a third-order RK time integration. On the other hand, the differences in the dispersive error between the two schemes is less significant, with the staggering scheme having a time shift around 2% of the initial wave period, compared to 4.7% for the HLLC scheme. It is worth noting that the values of both the diffusive and dispersive errors will vary depending on the wave parameters, along with the grid size.

For the following numerical tests, we use a water depth of $d = 1$ m and fixed dispersion coefficients, $\mu_1 = 0.00155$ and $\mu_2 = 0.0792$, corresponding to the limits of the intermediate water regime, which also relate $\lambda_1 = 25d$ and $\lambda_2 = 2d$. A fixed steepness coefficient of $\epsilon = 4 \times 10^{-5}$ is chosen to ensure that the wave remains within the linear
<figure><img src="image_7.png" /><figcaption>Fig. 4. Computed numerical diffusive and dispersive errors for a monochro- matic wave propagation after 20 cycles, for different water depths, while keeping the wave nonlinear and dispersive properties constant $\epsilon_0 = 4 \times 10^{-5}$ and $\mu_0 = 0.02$. Both the scaled diffusive and dispersive errors stay approximately constant for different problem scales (i.e., water depth).</figcaption></figure>

regime, allowing for comparison with Nwogu's exact linear equations. We run the same test using different grid resolutions: $\lambda_0/\Delta x = 20, 25, 30, 35, 40, \text{ and } 45$.

The computed dispersive and diffusive RMSEs for different grid sizes are shown in Fig. 5. It is important to note that the error shown in Fig. 5 is calculated after 20 cycles of wave propagation. Since both diffusive and dispersive errors are accumulative, large errors should be expected with a longer distance of propagation. All the numerical schemes show grid convergence for both the diffusive and dispersive errors. For the HLLC scheme, we note a significant reduction in the diffusive error and an improved convergence rate with higher-order reconstruction, while the dispersive error appears to be less sensitive to the reconstruction order. Overall, the conservative staggered scheme outperforms the HLLC scheme in terms of diffusive error, with the diffusive error being nearly an order of magnitude smaller than that of the HLLC scheme, despite the conservative scheme employing a discretization of only second-order accuracy. This can be attributed to the fact that the HLLC scheme and other Godunov-type schemes use an upwind approach for the discretization of the split part of the hydrostatic pressure term, while the staggered scheme applies a central difference method, which introduces less numerical diffusion. This can also explain the high dissipation error in the long wave regime ($\lambda_0 = 25d$) – especially for the HLLC scheme – where the hydrostatic pressure terms are of greater magnitude than in short and relatively non-hydrostatic waves. For the dispersive error, both the conservative scheme and the HLLC scheme exhibit comparable errors that decrease with grid size. The conservative

staggered scheme has a faster convergence rate for longer waves, while the HLLC scheme shows a faster convergence rate for shorter waves.

The dashed lines in Fig. 5 show the compound phase error, which include both the numerical error and dispersion error from Nwogu's equations. For very long wave, $\lambda_0 = 25d$, the dispersion error from the governing equations is negligible, causing the dashed and solid line to superpose. For a wave in intermediate water, with $\lambda_0 = 2d$, where the phase-speed from Nwogu's equations starts to overestimates that of Airy wave theory by 1%, the phase shift becomes noticeable — mostly due the dispersion error from the approximations in Nwogu's equations. Additionally, the sign of the numerical error affects the compound phase shift: Nwogu's equations, which have a positive dispersion error, tend to counteract the numerical error inherent to the conservative staggered scheme, while for the HLLC scheme the errors add up.

### 3.2. Spectral wave propagation

We further examine the impact of numerical errors, particularly diffusive errors, on the propagation of a wave spectrum, i.e. on an ensemble of various monochromatic waves with different amplitudes, frequencies, and phases. While for monochromatic waves only a single wave is affected by the diffusion error, which, theoretically, would allow for a careful choice of the grid size to minimize it, in real-world applications where waves of multiple frequencies are present, it is much more difficult, as numerical dissipation varies across different frequencies.

For the numerical test, we generate an irregular wave field based on the JONSWAP empirical spectrum with a significant wave height of $H_s = 1$ m and a peak period of $T_p = 15$ s in 50 m water depth, resulting in a peak wavelength of approximately $\lambda_p \approx 283$ m. The wave field is composed of 659 waves with a range of wavelengths between 100 m and 640 m. Each wave has a random phase between 0 and $2\pi$. For 2-hours simulation, the corresponding $\epsilon$ values range from $6 \times 10^{-3}$ to $8 \times 10^{-2}$, and $\mu$ values range from $2 \times 10^{-9}$ to $3 \times 10^{-5}$. These values indicate that the individual waves remain within the linear regime, as shown in Le Méhauté (1976). This ensures that nonlinear wave-wave interactions are negligible. Consequently, the wave spectrum is expected to conserve both its energy and shape while propagating over a flat bottom without significant exchanges among frequencies.

In a 15 km long flat channel, we record the free-surface at 10 km from the wavemaker. We run the test using the second-order Conservative Staggered scheme (both in space and time) and the HLLC scheme using a third-order time integration in combination with a fifth-order reconstruction and different grid sizes: $\Delta x = 2$ m, 5 m, 10 m, and 15 m. For each computation the Welch method is used to determine the power spectral density from a 6200-second free surface time-series, sampled at intervals of $\delta t = 0.5$ s.

Fig. 6 presents the computed spectral density distributions for two different schemes: the Conservative Staggered scheme with second-order integration in both time and space, and the HLLC scheme with third-order time integration combined with fifth-order flux reconstruction. Each plot represents a different grid size. For a grid size of $\Delta x = 2$ m, which corresponds to 141 grid points per peak wavelength, both schemes conserve the total wave energy and spectral shape after 10 km distance of propagation. The two lines representing the different numerical schemes, (blue for the conservative staggered and red for HLLC) are overlapping, which show both scheme converge to the exact solution when the grid is refined. It is important to note that some of the observed differences in the spectral distribution, such as the dip around the 0.09 Hz frequency, can be attributed to the spectral density estimation method, and not to the numerical solution. For a grid size of $\Delta x = 5$ m, which corresponds to 56 grid points per peak wavelength, we start to notice some differences between the two schemes, where the HLLC scheme start to have some dissipation of the energy in the high frequency tail (i.e., $f > 0.1$). This loss of energy in the high-frequency band becomes more noticeable as the grid size increases,
<figure><img src="image_8.png" /><figcaption>Fig. 5. Computed diffusive and dispersive errors after 20 cycles of wave propagation for different grid sizes and numerical schemes. Left: numerical error in a linear wave regime in intermediate water (i.e., $λ_0 = 2d$, corresponding to $kd = π$); Right: numerical error in a linear wave regime in shallow water (i.e., $λ_0 = 25d$, corresponding to $kd ≈ 0.25$). The dashed line represents the total error compared to Airy wave theory, this error combines both the dispersion numerical error from the schemes and the dispersion error from the Nwogu's equations. For $λ_0 = 25d$, the solid and dashed lines overlap, as the dispersion error from Nwogu's equations is negligible.</figcaption></figure>

with energy corresponding to frequencies associated with 15 grid points or fewer being completely absent from the spectrum. In contrast, the conservative staggered scheme shows less sensitivity to the grid size, with the overall wave energy reduced by only 0.5%, 1%, 3%, and 6% for $Δx = 2$ m, $Δx = 5$ m, $Δx = 10$ m, and $Δx = 15$ m, respectively. The modeled shorter waves are especially affected by numerical diffusion for two compound reasons. First, for a given grid size, shorter waves are represented by fewer grid points, as seen in Fig. 6. Second, over the same distance of propagation, shorter waves complete more cycles than longer waves. Since shorter waves are disproportionately impacted by the numerical dissipation, a more diffusive scheme, such as the HLLC scheme, results in a greater reduction in wave energy at higher frequency, thus distorting the shape of the spectrum. This can lead to changes in the minimum return group period ($T_{Rmin}$), potentially increasing it as high-frequency components diminish. Such a shift implies that energy-carrying wave groups arrive less frequently, which could alter the regularity and timing of wave events reaching the shore.

The execution time for computing the 2-hour spectral wave propagation with each numerical scheme is shown in Table 1. The numerical solution was implemented in C++ and parallelized using OpenMP (see the shared repository WavePropError (Mihami, 2025)). Model runs were performed on a 13th Gen Intel(R) Core(TM) i9-13900H processor with

20 cores. Both the Conservative Staggered scheme and HLLC scheme solutions share the same implementation for the wave-maker, sponge layer, and the discretization of the dispersion terms. Therefore, any differences in computational efficiency between the two schemes are primarily due to the numerical solution of the hydrostatic Shallow Water terms. For the same order of discretization (i.e., 2nd-order in space and time), the Conservative Staggered scheme is around 1.7 faster than the HLLC scheme for the same run, this is mainly due to its lean upwind flux computations, which avoid the approximate Riemann solution. Moreover, combining staggering of the variables in time with the Leapfrog method for the hydrostatic pressure term provides greater efficiency compared to the second-order Runge-Kutta method used in the HLLC scheme. As shown in the previous tests with monochromatic wave propagation, the HLLC scheme with only second-order discretization is highly diffusive, causing the wave amplitude to vanish completely after just a few cycles. It is important to note that the spectral results shown in Fig. 6, are computed using a 5th-order flux reconstruction combined with a 3rd-order time integration, which requires additional computational cost. For a converged numerical solution (i.e., $Δx = 2$ m), the model run using the high-order HLLC scheme takes approximately 480 s, nearly twice as long as the second-order HLLC scheme and over three times longer than the Conservative
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Spectral density distribution of the JONSWAP spectrum ($H_s = 1$ m and $T_p = 15$ s) in 50 m water depth after 10 km of propagation, computed using a second-order Conservative Staggered scheme (both spatially and temporally), and the HLLC scheme with third-order time integration and fifth-order reconstruction, for different grid resolutions. The x-axis represents the wave frequency and corresponding grid count per wavelength.</figcaption></figure>

<caption>**Table 1**</caption>

Execution times in seconds for a 2-hour spectral wave propagation in 15 km long channel for different numerical implementations: Conservative Staggered Scheme (Present scheme) with 2nd-order time and space discretization, HLLC scheme with 2nd-order time and space discretization, and HLLC scheme with 3rd-order time integration and 5th-order flux reconstruction. Simulations were performed on a 13th Gen Intel(R) Core(TM) i9-13900H processor with 20 cores.

<table><thead><tr><th rowspan="2">Grid size</th><th colspan="3">Execution time in seconds</th></tr><tr><th>Present scheme<br>2nd-order<br>in space/time</th><th>HLLC scheme<br>2nd-order<br>in space/time</th><th>HLLC scheme<br>5th-order flux recon.<br>3rd-order time</th></tr></thead><tbody><tr><td>2 m</td><td>142</td><td>248</td><td>479</td></tr><tr><td>5 m</td><td>24</td><td>50</td><td>76</td></tr><tr><td>10 m</td><td>12</td><td>20</td><td>27</td></tr><tr><td>15 m</td><td>5</td><td>9</td><td>10</td></tr></tbody></table>

Staggered scheme. This highlights the numerical efficiency of the conservative staggered grid, which requires less computational time to achieve the same level of accuracy as the HLLC scheme. For coarser grids, the differences in computational speed between the different numerical methods are much smaller. This is largely because the time for the actual computation is rather short (under 1 min) in relation to the overhead from reading and saving data. It is worth noting that while the absolute time difference between the two numerical schemes may appear minor for this test (order of a few minutes), in

two-dimensional large-scale scenarios, these differences can scale up to hours of reduction in computation time.

The key findings from this analysis can be summarized as follows:

* **Numerical errors in linear wave propagation:** Both the Conservative Staggered scheme and the HLLC scheme introduce numerical errors that affect wave phase and amplitude, which can be mitigated by mesh refinement.

* **Differences in Numerical Schemes:** While both schemes use central finite differences for dispersion terms, the Conservative
Staggered scheme applies central differences for the hydrostatic pressure (i.e., the free surface gradient $gH\eta_x$), whereas the HLLC scheme splits this term into flux and bottom gradient terms, using upwind methods that introduce numerical diffusion, which in turn can be reduced with higher-order flux reconstruction.

*   **Error Characteristics:** The Conservative Staggered scheme generally tends to underestimate wave celerity and wavelength (negative dispersion error), while the HLLC scheme overestimates them (positive dispersion error). The dispersion error magnitude is comparable but varies with wave frequency and water depth.

*   **Sensitivity and Accuracy:** The Conservative Staggered scheme is less sensitive to diffusion errors and requires fewer grid points to maintain wave energy and frequency distribution compared to the higher-order HLLC scheme.

*   **Computational Efficiency:** The lean numerical structure of the Conservative Staggered scheme results in greater computational efficiency and faster performance than the Riemann-based HLLC scheme.

## 4. Sensitivity of the staggered scheme to wave-breaking

Excessive numerical diffusion is generally an undesirable aspect of numerical methods as it can reduce the accuracy of the solution by smoothing out critical wave features. However, a controllable amount of dissipation is often necessary to ensure stability and convergence — especially in cases that involve wave breaking. Since the standard Boussinesq-type formulation is derived for non-breaking waves, a wave-breaking closure needs to be implemented to account for the effects of dissipation and to prevent from singularities as the mesh is refined. From a large-scale perspective, the main consequence of wave breaking is strong energy dissipation. One way to mimic these effects is by introducing a diffusion term into the momentum equations, where eddy viscosity serves as a wave-breaking proxy that quantifies the effective viscosity caused by turbulent fluctuations in the flow (see Zelt (1991)). Another more basic strategy, is the “hybrid approach” (Kazolea and Ricchiuto, 2018), which involves deactivation of the dispersion terms in the governing equations. This technique subsequently reduces the solution locally and momentarily to that of the hyperbolic shallow water equations and the approximation of wave breaking as a propagating bore or discontinuity. A variety of shock-capturing numerical schemes can support these properties of the shallow water equations (Toro, 2001), which has been a key argument for the shift from finite difference (FD) methods to hybrid finite volume-finite difference (FV-FD) schemes in the discretization of Boussinesq-type equations (e.g., Roeber et al. (2010), Shi et al. (2012)). However, the accuracy and stability of the approach depend on numerical dissipation to bypass the inherent incompatibility between the dispersive and locally hydrostatic solutions. This means that the solution can be very sensitive to the grid size as shown in Kazolea and Ricchiuto (2018). In contrast, the eddy viscosity approach has proven to be less sensitive to the grid size (see Kazolea and Ricchiuto (2018)) — especially when it is based on a spatially and temporally varying TKE approximation. This makes it an attractive choice for low-dissipative numerical methods such as the conservative staggered scheme.

In this second part, we extend the validation of the conservative staggered scheme for cases that involve wave breaking. We assess the accuracy and stability of the TKE-based and hybrid wave-breaking approaches when combined with the conservative staggered scheme and analyze their performance in both solitary wave and spectral wave-breaking scenarios.

### 4.1. Wave breaking of a solitary wave

In this initial test, we focus on validating the numerical scheme for modeling wave-breaking processes of a solitary wave. For that, we use the laboratory data from the HReef experiment conducted by Roeber

et al. (2010), Roeber and Cheung (2012) at the O.H. Hinsdale Wave Research Laboratory at Oregon State University. This benchmark test has been widely used to validate Boussinesq-type models, as it captures several key wave processes, such as shoaling, bore formation, wave reflection, and swashing (see Roeber and Cheung (2012), Filippini et al. (2016), Kazolea and Ricchiuto (2018), Tonelli and Petti (2011)). The flume used in the experiment has an effective length of 83.7 m, with a 1:12 slope starting at 25.9 m from the wavemaker. At a water level of 2.5 m, the lagoon area has a shallow depth of 0.136 m, while the reef crest remains exposed by 0.065 m. Additional details of the experiment setup can be found in Roeber et al. (2010). The benchmark test is computed in a 1D computational domain using a Courant number $Cr = 0.5$ and a Manning coefficient of $n = 0.013 \text{ s/m}^{1/3}$ to describe the smooth concrete and plywood surfaces of the flume. A 0.75 m solitary wave is generated at the left boundary of the domain by imposing the following boundary condition:

$$ \eta(x = 0, t) = a_0 \operatorname{sech}^2 \left[ \sqrt{\frac{3a_0}{4h_0^3}} \left( (t_0 - t) \sqrt{g(a+h_0)} \right) \right] \\ u(x = 0, t) = \frac{\eta(x = 0, t)}{h_0 + \eta(x = 0, t)} \sqrt{g(a_0 + h_0)}. $$

where $h_0 = 2.5$ m is the water depth at the offshore boundary, $a_0 = 0.7$ m is the amplitude of the solitary wave, and $t_0 = 41.52$ s is the time shift used to align the solitary wave crest with the laboratory data. Once the solitary wave is generated, the boundary condition is changed to a reflective one.

The numerical solution is computed with a grid size $\Delta x = 0.05$ using the conservative staggered scheme with TKE-based wave breaking closure as described in Section 2. The resulting free surface and TKE values are shown in Fig. 7. As shown in the figure, the solitary wave undergoes several key processes:

*   **Energetic wave breaking** as the solitary wave shoals and interacts with the exposed reef. The wave forms a large air cavity over the dry reef crest, with the front face impinging into the lagoon.

*   **Propagation and reflection of the bore** across the lagoon.

*   **Generation of an undular bore** (or a train of dispersive waves) as the reflected bore enters the deeper water offshore of the reef crest.

The numerical solution is able to correctly account for the plunging breaker effects using the eddy viscosity approach. At the wavefront of the breaking wave, the computed TKE values increase sharply, creating a localized diffusion zone that mimics physical turbulent dissipation. The TKE is then generated and advected, following the hydraulic jumps as it moves and reflects within the lagoon. As the reflected wave overtops the reef crest and propagates from shallow to deeper water, the effect of the eddy viscosity is reduced, allowing the transition of the bore into a series of dispersive waves. It is worth noting that the computation of the moving boundaries with the present scheme requires no additional treatment, yet it provides good stability and robustness as shown in this test and also in a more detailed run-up study in Mihami et al. (2022).

To gain more insight into the effectiveness of the TKE-based wave-breaking closure in the context of the conservative staggered scheme, the test is recomputed using the hybrid wave-breaking approach commonly employed in many Boussinesq-type models. In this case, the dispersion terms are locally and temporarily deactivated so that the numerical solution represents only the SWE in particular cells where the following breaking criterion is met:

$$ \frac{\eta}{h} > C_b, \quad \text{with} \quad C_b = 0.78. \qquad (37) $$

Although simplified, this breaking criterion was shown to be the most robust when used with the hybrid approach, as reported in Shi et al. (2012), Tonelli and Petti (2009).
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Snapshots of free surface profiles for the propagation of solitary wave with $A/h = 0.3$ over 1:12 slope and exposed reef crest. Blue lines and red lines denote the water free surface and the TKE. The TKE vertical axis is zoomed in to highlight the shape of the TKE solution. During the computation, the TKE quantity reaches a maximum value of $33 \text{ m}^2 \text{s}^{-2}$ as the wave begins to break around 31 s, when the wavefront reaches $x = 50 \text{ m}$.</figcaption></figure>

**Fig. 8** compares the time series from the computed and recorded surface elevation. While the TKE-based eddy viscosity approach implemented in the model allows for an accurate and stable description of the surging and plunging waves over the reef with a relatively fine grid $\Delta x = 0.05 \text{ m}$, the numerical solution with the hybrid breaking closure becomes unstable with $\Delta x = 0.05 \text{ m}$. The instabilities begin at $t = 55 \text{ s}$, which is equivalent to $t\sqrt{g/h} = 109$. This occurs as the reflected bore propagates from the reef to deeper water, releasing super-harmonics.

The sudden transition from the hydrostatic to the dispersive solution is what causes the numerical solution to become unstable. However, the solution remains stable when the hybrid approach is applied with a coarser grid size of $\Delta x = 0.1 \text{ m}$, as shown by the red lines in Fig. 8, this is mainly due to the increased numerical diffusion introduced by the coarser grid, which suppresses the growth of instabilities. A minor spurious oscillation is observed at gauge $x = 54.4 \text{ m}$ around
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Time series of the free surface at different wave gauges for the propagation of solitary wave with $A/h = 0.3$ over 1:12 slope and exposed reef crest. The black circles represent the measured data. The blue solid line corresponds to the computed results with a TKE-based eddy viscosity approach with a grid size $\Delta x = 0.05$ m. The red line shows the solution obtained with a hybrid wave-breaking closure with a grid size $\Delta x = 0.1$ m.</figcaption></figure>

$t\sqrt{g/h} = 120$; however, this oscillation is quickly damped and does not significantly affect the numerical solution.

With the TKE approach, the additional diffusion term, which is dependent on the computation of TKE, effectively prevents from potential instabilities, which often arise from depth-integrated dispersive solutions in combination with sharp flow transitions. The dissipation has only a local impact on the numerical results. This means that the dispersion properties of the equations are not altered, and subsequent wave transformation processes, such as the release of high-order harmonics, are well accounted for.

## 4.2. Breaking of irregular waves

In this section, we extend the validation of the conservative staggered scheme to include wave breaking in spectral waves. Unlike

solitary wave tests, spectral waves consist of multiple frequencies, resulting in less localized breaking that occurs intermittently across the surf zone. For this validation, we use the experimental data from the WEN_004 test presented in Suzuki et al. (2017). This test was conducted in a wave flume with a length of approximately 70 m, a width of 4 m, and a uniform water depth of 0.942 m in the offshore region. The flume included a sloping beach with a 1:35 gradient, designed to induce wave breaking as the waves propagate toward the shore. A piston-type wave generator produced irregular waves based on a JONSWAP spectrum with a peak period of 11.75 s and a significant wave height of 0.093 m. Detailed measurements of the free surface elevation were recorded using an array of 12 wave gauges distributed along the flume. Although the test was conducted to study overtopping, the data can be used to validate the numerical method in simulating wave
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Numerical model layout showing the location of the gauges along with the position and width of the wavemaker and sponge layer for wave generation.</figcaption></figure>

shoaling across the full wave spectrum, as well as the breaking process and the subsequent transfer and dissipation of wave energy.

To accurately reproduce the benchmark test with the correct wave input, including phase information, the free surface elevation recorded by the first wave gauge (WG1) in the experiment is used as wave input. This is implemented through internal wave generation, as described by Wei et al. (1999), with a sponge layer condition applied on the left side of the domain to absorb incoming reflected waves. This method has been chosen over the standard boundary condition input because it addresses the issue of missing velocity data from the laboratory measurement, which is required for the model.

To accurately generate the input wave signal using the internal wavemaker, the time series recorded at the first wave gauge in the flume is analyzed in the frequency domain. This is done by applying a Fast Fourier Transform (FFT), which decomposes the free surface elevation signal into its constituent frequency components, along with their corresponding amplitudes and phase angles. Only frequency components corresponding to nondimensional wavenumbers in the range $0.2\pi \le kh \le 2\pi$ are retained. The high-frequency cutoff ensures that all wave components fall within the validity range of Nwogu's dispersion relation, thus preserving the accuracy of the wave model. The low-frequency cutoff excludes long-wave components that contribute less than 1% of the total energy, thereby avoiding the need for an excessively wide wave-maker or sponge layer to accommodate low-frequency, high-wavelength waves.

After applying these cutoffs, the amplitude spectrum is uniformly rescaled so that the significant wave height ($H_s$) of the modified signal matches the original. This means that the energy removed at low and high frequencies is redistributed evenly across the retained frequency components, ensuring that the total wave energy remains consistent.

To validate the numerical method, we first run the test with a grid size $\Delta x = 0.05$ m using the TKE eddy viscosity approach. A Manning's coefficient of $n = 0.012 \text{ s/m}^{1/3}$ was selected to represent the smooth wood used in the experiment. The wave maker is positioned at $x = 6$ m with a width of $\lambda_{max}/2$, where $\lambda_{max} = 9.42$ m is the longest wave of the reconstructed spectrum. The sponge layer width is set to $\lambda_{max}/4$ to ensure optimal absorption (see Fig. 9). It is important to note that the phases for the wavemaker were calculated from the measured time series at the first gauge, WG1 and it is therefore expected that the time series closely matches the measured signal at WG1.

Fig. 10 shows the time series of free-surface elevation recorded at different wave gauges: WG3, WG5, WG7, WG8, WG10, and WG12. The results indicate a close agreement between the numerical solution and the experimental data, where the overall wave propagation, shoaling, breaking, and overtopping/reflection processes are matched. A small time shift of approximately 0.5 s can be observed at the gauges closer to the beach, starting from WG8. Though this shift is consistent with the numerical dispersion error inherent in the conservative staggered scheme, it appears to be independent of the grid resolution. Even with

a finer grid resolution of $\Delta x = 0.02$ m, the phase shift remains. This suggests that the shift is connected to the wave-breaking closure and the threshold chosen for the production of the TKE, which is still able to describe the overall wave-breaking process despite the limitations.

In this second part, we compare the model results for different wave-breaking closure methods: (1) using the TKE-based eddy viscosity approach, (2) applying the hybrid approach, which relies on deactivating the dispersion terms and switching from a dispersive to hydrostatic solution, and (3) without any explicit wave-breaking treatment, instead relying on numerical dissipation. For each method, we compare the results for two grid sizes: $\Delta x = 0.12$ m and $\Delta x = 0.06$ m, the coarser grid to highlights the effects of numerical dissipation on wave breaking and the finer grid to assess the converged solution. The benchmark test is run for 2400 s, with numerical results saved at a sampling interval of $dt = 0.05$ s for each gauge. It is worth noting that the solution using the hybrid approach with any grid size smaller than $\Delta x = 0.06$ m resulted in instabilities, preventing the simulation from running for the full 2400 s. These instabilities arise at the interface where the transition occurs from the dispersive to the hydrostatic solution (see Appendix A.4).

Fig. 11 presents the spectral distribution computed with both the coarse and fine grids at four wave gauges. For the coarser grid ($\Delta x = 0.12$ m), the differences between the different wave-breaking methods are generally small and show good agreement with the laboratory data. This highlights the ability of numerical dissipation to effectively capture some of the energy transformation associated with wave breaking. However small discrepancies can be still observed at WG8, particularly at the swell peak $f = 0.4$ Hz, where the no-breaking model over-predicts energy, while the hybrid and TKE-based approaches show better agreement. As the grid size is refined, the differences between the wave-breaking methods become more pronounced. In the absence of a wave-breaking closure, the model overestimates the wave energy across all frequencies in the surface zone (WG8 and WG12). The TKE-based approach provides the best overall agreement with the experimental data, particularly at higher frequencies, where it effectively dissipates excessive energy. The hybrid approach, while still improving the high-frequency energy prediction compared to the no-breaking model, tends to overestimate energy in the infragravity band for all the gauges. We attribute these artificial IG waves to the abrupt transition from dispersive to hydrostatic regime in the hybrid approach, where dispersion is switched off locally and instantaneously once the breaking criterion is satisfied. This sudden change can introduce spurious long-period components, particularly on refined grids where the transition is more pronounced (see Kazolea and Ricchiuto (2018)). In contrast, such artifacts are not observed in the TKE-based closure, which introduces the dissipation gradually through a dynamically adjusted eddy viscosity term. This suggests that while numerical dissipation can play a role in capturing wave-breaking effects on coarser grids, explicit turbulence-based closures are necessary to achieve accurate energy dissipation at finer resolutions.
<figure><img src="image_13.png" /><figcaption>Fig. 10. Time series of the free surface at different wave gauges. Black lines denote the measured data, while red lines denote the numerical model results, obtained with a grid size of Δx = 0.05 m and using the TKE-based wave-breaking approach.</figcaption></figure>

Fig. 12 extends the analysis by investigating the sensitivity of wave-breaking methods to grid resolution. Unlike the previous test, which focused on two specific grid sizes (Δx = 0.12 m and Δx = 0.06 m), this analysis systematically varies the grid size from Δx = 0.15 m to Δx = 0.01 m in increments of 0.01 m. This allows us to assess how each wave-breaking method responds to grid refinement and whether the results converge to a stable solution. The computed significant wave heights (H_s) are compared against experimental measurements at each gauge, with marker intensity representing grid resolution — darker markers for coarser grids and lighter markers for finer ones. The results highlight the natural tendency of the original Boussinesq-type equations (without any wave-breaking closure) to overestimate wave energy in the surf zone when the grid is refined. There is also a noticeable overestimation of the wave height during shoaling, particularly at gauge WG7. This is consistent with the finding in Filippini et al. (2015), where it was shown that weakly nonlinear equations, such as Nwogu's equations, tend to overestimate the peaks of the shoaling waves compared to fully nonlinear equations. The hybrid breaking approach shows a high sensitivity to grid size, with instabilities occurring for fine grids with a resolution of Δx = 0.05 m or smaller. The TKE wave-breaking method, on the other hand remains stable for all grid sizes up to Δx = 0.02 m. To further quantify the sensitivity of each wave-breaking method to grid resolution, we define a metric δ_H_s as the mean standard deviation of the significant wave height across all tested grid sizes at each gauge:

$$ \delta_{H_s} = \frac{1}{M} \sum_{j=1}^{M} \sqrt{\frac{1}{N} \sum_{i=1}^{N} (H_{s,i}^{(j)} - \bar{H}_{s}^{(j)})^2} $$

where H_{s,i}^{(j)} is the computed significant wave height at gauge j for grid size i, H_s^{(j)} is its mean over all grids, N is the number of grid sizes, and M is the number of gauges. This measure reflects the variability of every breaking method with respect to grid refinement.

The results show that the original Boussinesq-type equations (without wave breaking) exhibit the largest variability (δ_H_s = 0.006 49 m), reflecting their strong dependence on grid size and their tendency to overestimate wave energy in the surf zone. The hybrid breaking approach shows moderate variability (δ_H_s = 0.003 34 m), particularly for finer grids where the over-shoaling from Nwogu's equations becomes more pronounced and the transition from dispersive to hydrostatic solutions becomes increasingly difficult. In contrast, the TKE wave-breaking method achieves the lowest variability (δ_H_s = 0.001 77 m), confirming its robustness and stability across different grid resolutions, and providing a good match with experimental data. This makes TKE a promising approach to address the over-shoaling issues associated with weakly nonlinear Boussinesq equations, offering a more robust solution for simulations with refined grids.

### 4.3. Wave breaking over 3D reef

Having validated the conservative staggered scheme for both solitary and spectral wave breaking for two-dimensional tests, we now extend our validation to a more complex bathymetric setting: wave breaking over a three-dimensional reef, using the experimental data from Swigler (2010), Lynett et al. (2010). Unlike previous cases, where wave evolution was primarily governed by shoaling and wave-breaking
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Power spectral density of the measured and computed free surfaces at four wave gauges. Results are shown for two grid resolutions: a coarse resolution ($\Delta x$ = 0.12 m, left column) and a fine resolution ($\Delta x$ = 0.06 m, right column). Three numerical approaches for wave breaking are compared: (1) no breaking implementation (red), (2) hybrid breaking approach (green), and (3) TKE-based breaking approach (blue).</figcaption></figure>

over relatively simple bathymetries, this test involves additional 3D wave processes such as wave refraction, diffraction, and overtopping. Moreover, the test serves as a useful benchmark for assessing the robustness of the wetting-drying algorithm. This validation test was specifically chosen due to its complexity, since it was mentioned to present stability challenges in traditional collocated finite-difference methods, as highlighted in Shi et al. (2012).

Fig. 13 shows the laboratory basin obtained from a LiDAR scan, along with the layout of the gauges, which include resistance-type wire wave gauges (positions 1 to 9) and ADVs (positions 2, 3, and 10 only). The bathymetry consists of a planar beach with a 1:30 slope, connected to a triangular-shaped shelf and a conical island located on the shelf edge. The numerical computation is set up with a grid size of $\Delta x = \Delta y = 0.1$ m, a Courant number of $Cr = 0.5$, and a Manning friction coefficient $n = 0.014 \text{ m/s}^{1/3}$. This setup is consistent with those used in previously published studies (e.g., Roeber and Cheung (2012), Shi et al. (2012)). The incident solitary wave is generated at the boundary

with a height of 0.39 m in a 0.78 m water depth, resulting in strongly nonlinear conditions with $A/h = 0.5$. For wave breaking, we employ the TKE-based eddy viscosity approach, with the Froude number at the free surface serving as criterion for the onset of TKE production. The computed results of the free surface for the entire domain are shown at selected instances in Fig. 14 for illustration. Since more details about the wave processes can be found in Lynett et al. (2010), the following description mainly focuses on the TKE-related processes.

Upon wave breaking, TKE is generated in both the x- and y-directions and is advected with the breaking wave fronts to either side of the apex. This process ensures that the dissipation related to wave breaking is accurately represented while maintaining the dispersive nature of the solution. At $t = 5.4$ s, the solitary wave reaches the shelf apex, breaking as a strong plunging breaker due to the steep slope. The computed TKE peaks at $0.527 \text{ m}^2/\text{s}^2$ along the centerline just onshore of the apex. The breaking bore then expands from the centerline to the basin side wall, where the peak TKE values of the simulation reach
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Significant wave height (H_s) at different wave gauges. Computed results are presented for three different wave-breaking models: (1) without any breaking model (red), (2) with the hybrid breaking model (green), and (3) with the TKE-based breaking model (blue). The results also cover a range of grid sizes from Δx = 0.01 m to Δx = 0.15 m, with the grid size indicated by the color shading of the marker: darker shades represent coarser grids, and lighter shades correspond to finer grids. The measurements (black line) are shown alongside the computed for each plot.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Breaking solitary wave over 3D reef. Layout with locations for resistance-type wire wave gauges and ADVs (positions 2, 3, and 10 only).
*Source:* Adapted from Roeber and Cheung (2012).</figcaption></figure>

15.5 m²/s² at t = 7.9 s, driven by intense shoaling over the gentler parts of the slopes near the basin side walls. A secondary bore develops at the shelf and moves offshore as the flow relaxes. TKE values on the lee side of the cone, observed from t = 7.80 s to t = 12 s, highlight the offshore bore, which is consistent with the experiment findings (Lynett et al., 2010). As the wave reaches the dry area, small vortices are generated as the main bore is reflected off the beach, with TKE values emerging in response to the energy associated with vortex formation.

**Fig. 15** compares the computed and recorded surface elevations over the full trial period of 90 s. The model reproduces the measured surface elevations in front of and behind the cone, where wave breaking and overtopping happen. As shown in Fig. 16, the numerical solution also matches the cross- and longshore-components of the velocity reasonably well. It should be noted that the velocities in the physical model were measured near the free surface whereas the numerical solution returns the velocity close to mid-depth. The overall agreement between the computed and recorded data demonstrates the validity of the proposed formulation in handling multiple hydraulic processes, transitions from hydrostatic to dispersion-dominated flows, as well as a variety of wave-breaking scenarios in the two-dimensional horizontal plane with a moving boundary.

## 5. Summary and discussion

This paper presents and validates a conservative staggered-grid scheme for the numerical solution of weakly dispersive Boussinesq-type equations, with particular emphasis on the formulation proposed by Nwogu (1993). The proposed approach re-examines the core discretization strategy – specifically the use of a momentum-conserving staggered grid layout – to improve both numerical accuracy and stability. This reformulation thereby questions the prevailing reliance on hybrid finite volume–finite difference (FV-FD) schemes built around Riemann-based solvers. The results demonstrate that a carefully designed conservative staggered scheme can significantly reduce numerical diffusion and computational overhead, while preserving the accuracy and stability of wave transformation and breaking processes.
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Evolution of the computed free surface (shown in the upper 3D plots) and the turbulent kinetic energy (TKE) values (shown in the lower surface plot) at various time steps throughout the simulation: - (a) *t* = 5.4 s: Solitary wave begins breaks at the apex of the reef flat; (b) *t* = 7.9 s: Wave overtops the cone and continue propagating as a breaking bore; peak TKE values near the lateral boundaries of the domain; (c) *t* = 9.0 s: Wave continues breaking as it approaches dry beach with high TKE values observed around the cone as the wave wraps around it; (d) *t* = 19.8 s: Run-up on the dry concrete beach with sloshing and small vortices moving around the cone.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Free surface elevation at the gauges. The blue line denotes the computed result and the black circles represent the laboratory data.</figcaption></figure>

The numerical tests highlight several key advantages of the proposed scheme. First and most notably, in the context of long-distance

wave propagation, the scheme exhibits less sensitivity to grid spacing, which results in substantially lower numerical dissipation compared to
<figure><img src="image_19.png" /><figcaption>Fig. 16. Horizontal flow velocities at the gauges. The blue line denotes the computed result (corresponding to mid-depth), the black circles represent the laboratory data (measured near the free surface).</figcaption></figure>

conventional HLLC-based solutions. The key advantage is reflected in the fact that significantly fewer grid cells per wavelength are needed to achieve nearly diffusion-free wave propagation. This is particularly evident in both monochromatic and spectral wave propagations, where the conservative staggered method preserves wave amplitude and spectral shape with only a fraction of the grid resolution typically required by Godunov-type solvers. In fact, the diffusion error was found to be one order of magnitude smaller than that of the HLLC scheme, despite operating at only second-order accuracy in both space and time. It is important to note, however, that both the conservative staggered and the high-order HLLC schemes converge to nearly identical solutions when sufficiently high resolution is used. This reinforces that the HLLC approach is not inherently incorrect or unsuitable for nearshore wave modeling, but rather less practical due to its high demands on grid resolution. In many real-world coastal modeling applications, the number of grid cells per wavelength falls well below the threshold required for accurate HLLC performance, simply because the associated computational cost becomes impractical in large-scale or time-sensitive simulations. Second, while the dispersive error introduced by the proposed scheme is comparable to that of the HLLC scheme, its behavior offers a practical advantage. Specifically, the conservative staggered scheme introduces a slight negative phase error, which partially offsets the intrinsic positive dispersion error associated with the (Nwogu, 1993) formulation. Though not guaranteed for all Boussinesq-type equations, this fortuitous cancellation can result in improved overall phase accuracy, particularly in the intermediate water depth regime where dispersive effects are most pronounced. It is worth to highlight that the proposed second-order conservative staggered scheme was benchmarked against a high-order HLLC solver employing fifth-order flux reconstruction and third-order Runge-Kutta integration. Despite its lower formal order, the staggered scheme achieves comparable or superior accuracy with a fraction of the computational cost. Its lean formulation – avoiding Riemann solvers and higher-order reconstructions – offers clear advantages for computational efficiency. However, it should be acknowledged that the staggered-grid implementation requires additional coding effort due to the staggering of the variables (e.g., interface vs. cell center), which involves careful handling of spatial indices and derivative computations

to ensure consistency. Boundary conditions are also less straightforward to implement compared to the HLLC framework, which uses a collocated arrangement.

Building on these results, the scheme's performance for modeling wave breaking and other nonlinear processes has been evaluated. Compared to earlier non-conservative staggered and collocated finite-difference (FD) approximations, the scheme demonstrates substantially improved stability, particularly in challenging scenarios that involve wave breaking, overtopping, and complex wave transformation. Nevertheless, the results confirm that the scheme performs best when combined with a turbulent kinetic energy (TKE)-based eddy viscosity model, rather than relying on standard hybrid wave-breaking approaches that deactivate dispersion terms. In contrast, the HLLC-based solver, while more diffusive, is less sensitive to the choice of the wave-breaking model. The higher background diffusivity tends to smooth out steep gradients, improving numerical robustness. However, this behavior is not advantageous, as it effectively relies on non-controllable numerical diffusion to stabilize breaking waves — a strategy that is case-dependent and difficult to generalize across different wave conditions.

Over a range of standard benchmark tests that include solitary wave breaking, irregular spectral waves, and three-dimensional reef cases — the TKE-based closure consistently provides smooth, localized, and physically consistent dissipation, especially under highly energetic breaking conditions. Its formulation, governed by the vertical gradients of horizontal velocities, enables a realistic and stable representation of surf zone dynamics without introducing spurious oscillations or compromising the dispersive properties of the solution. The results also highlight the TKE approach as a reliable solution to the over-shoaling issues commonly associated with weakly nonlinear Boussinesq-type models. Overall, the combination of the conservative staggered scheme with the TKE-based approach offers a robust, accurate, and scalable framework for high-fidelity coastal and nearshore wave modeling.

Despite these strengths, some limitations should be acknowledged

First, the turbulence model parameters, particularly those related to the eddy viscosity closure, remain empirically calibrated, and their generalization across a broad range of wave conditions requires further
study. This is especially true for the mixing length approximation, which has a significant influence on both the accuracy and the stability of the numerical solution. In this study, the mixing length was approximated as the still water depth, following recommendations from previous studies; however, a systematic investigation into the sensitivity of the model to different mixing length formulations could offer valuable insights. Additionally, the use of a two-equation turbulence closure model could potentially improve the consistency and physical realism of wave-breaking representations across diverse conditions, although this would inevitably introduce additional computational costs that must be carefully balanced.

Another limitation lies in the scope of the error analysis, which focused exclusively on linear wave propagation in order to enable direct comparison with the linearized analytical solutions of Nwogu (1993). While necessary for isolating numerical errors under controlled conditions, a more comprehensive investigation incorporating nonlinear wave propagation and higher-order interactions would be necessary to further validate the proposed scheme.

Finally, while the present work has relied on idealized analytical and experimental benchmark tests to validate the scheme's performance, the evaluation of the conservative staggered scheme in realistic coastal environments remains an important avenue for future research. Such applications would help quantify the practical advantages of the conservative staggered approach in operational settings, particularly in terms of accuracy, stability, and runtime performance.

In conclusion, the conservative staggered scheme presented in this study offers a compelling alternative to conventional FV-FD solvers for Boussinesq-type wave modeling. Its superior performance in minimizing numerical diffusion and reducing computational overhead – together with the incorporation of a robust, physically consistent wave-breaking closure – makes it a promising approach for improving the accuracy and efficiency of coastal and nearshore wave models.

### CRediT authorship contribution statement

**Fatima-Zahra Mihami:** Writing – original draft, Writing – review & editing, Visualization, Validation, Software (main), Methodology, Investigation, Formal analysis, Data curation, Conceptualization.

**Volker Roeber:** Writing – review & editing, Visualization, Validation, Software (support), Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

This research was carried out in the framework of the I-SITE program Energy & Environment Solutions (E2S). The authors acknowledge financial support from E2S-UPPA, the Communauté d'Agglomération Pays Basque (CAPB), and the Communauté Région Nouvelle Aquitaine (CRNA) for the chair HPC-Waves. This work also received additional support from the joint laboratory KOSTARISK, co-funded by E2S-UPPA, the AZTI Foundation, and the Rivages Pro Tech center of SUEZ. The authors further acknowledge funding from the European Union's H2020 programme under grant agreement No 883553. Dr. Volker Roeber gratefully acknowledges support from the European Union's Horizon Europe Marie Skłodowska-Curie grant programme under grant agreement No 101119437 and the Regions4Climate (R4C) – Building Resilient Communities EU project, grant agreement ID 101093873, funded under the HORIZON 2.5 Climate, Energy and Mobility programme (Call: HORIZON-MISS-2021-CLIMA-02-04).

At the University of Hawai'i, Dr. Fatima-Zahra Mihami was supported by the U.S. Office of Naval Research [Office of Naval Research] (ONR Award # N00014-23-1-2015); the United Nations Environment Programme/Green Climate Fund (UNEP/GCF Ref. # FP147); and the U.S. National Oceanic and Atmospheric Administration Award # NA21NOS0120091 to the Pacific Islands Ocean Observing System (PacIOOS). We would also like to thank Dr. Tomohiro Suzuki for sharing the laboratory data for the benchmarking test in chapter 4.2.

## Appendix

### A.1. Discretization of eddy viscosity terms

The eddy viscosity term is computed at the cell center $v_{i(i,j)}$ and wave breaking source term $[0, R_x, R_y]^T$ is approximated on a staggered grid as follows:

$$ 
\begin{aligned}
[R_x]_{i+\frac{1}{2},j} = & \frac{1}{\Delta x^2} \left[ v_{i(i+1,j)} H_{i+1,j} \left( u_{i+\frac{3}{2},j} - u_{i+\frac{1}{2},j} \right) \right. \\
& \left. - v_{i(i,j)} H_{i,j} \left( u_{i+\frac{1}{2},j} - u_{i-\frac{1}{2},j} \right) \right] \\
& + \frac{1}{2\Delta y^2} \left[ v_{i(i+\frac{1}{2},j+\frac{1}{2})} H_{i+\frac{1}{2},j+\frac{1}{2}} \left( u_{i+\frac{1}{2},j+1} - u_{i+\frac{1}{2},j} \right) \right. \\
& \left. - v_{i(i+\frac{1}{2},j-\frac{1}{2})} H_{i+\frac{1}{2},j-\frac{1}{2}} \left( u_{i+\frac{1}{2},j} - u_{i+\frac{1}{2},j-1} \right) \right] \\
& + \frac{1}{2\Delta x\Delta y} \left[ v_{i(i+\frac{1}{2},j+\frac{1}{2})} H_{i+\frac{1}{2},j+\frac{1}{2}} \left( v_{i+1,j+\frac{1}{2}} - v_{i,j+\frac{1}{2}} \right) \right. \\
& \left. - v_{i(i+\frac{1}{2},j-\frac{1}{2})} H_{i+\frac{1}{2},j-\frac{1}{2}} \left( v_{i+1,j-\frac{1}{2}} - v_{i,j-\frac{1}{2}} \right) \right]
\end{aligned}
$$

$$ 
\begin{aligned}
[R_y]_{i,j+\frac{1}{2}} = & \frac{1}{\Delta y^2} \left[ v_{i(i,j+1)} H_{i,j+1} \left( v_{i,j+\frac{3}{2}} - v_{i,j+\frac{1}{2}} \right) \right. \\
& \left. - v_{i(i,j)} H_{i,j} \left( v_{i,j+\frac{1}{2}} - v_{i,j-\frac{1}{2}} \right) \right] \\
& + \frac{1}{2\Delta x^2} \left[ v_{i(i+\frac{1}{2},j+\frac{1}{2})} H_{i+\frac{1}{2},j+\frac{1}{2}} \left( v_{i+1,j+\frac{1}{2}} - v_{i,j+\frac{1}{2}} \right) \right. \\
& \left. - v_{i(i-\frac{1}{2},j+\frac{1}{2})} H_{i-\frac{1}{2},j+\frac{1}{2}} \left( v_{i,j+\frac{1}{2}} - v_{i-1,j+\frac{1}{2}} \right) \right] \\
& + \frac{1}{2\Delta x\Delta y} \left[ v_{i(i+\frac{1}{2},j+\frac{1}{2})} H_{i+\frac{1}{2},j+\frac{1}{2}} \left( u_{i+\frac{1}{2},j+1} - u_{i+\frac{1}{2},j} \right) \right. \\
& \left. - v_{i(i-\frac{1}{2},j+\frac{1}{2})} H_{i-\frac{1}{2},j+\frac{1}{2}} \left( u_{i-\frac{1}{2},j+1} - u_{i-\frac{1}{2},j} \right) \right]
\end{aligned}
$$

### A.2. Solitary wave propagation

Solitary-wave propagation is a standard verification for Boussinesq-type models: with an exact initial profile, the wave should preserve shape, speed, and amplitude by nonlinear-dispersive balance. For Nwogu's equations, a physically admissible exact solitary wave is not available: the closed-form solution of Hamdi et al. (2005) requires $\sigma = -0.531$ and yields non-physical negative total depth, incompatible with conservative schemes. We therefore use the solitary wave from the approximate ODE reduction of Wei and Kirby (1995) to Nwogu's system (Eqs. (1)–(4)), which provides a practical, admissible benchmark; see Wei and Kirby (1995) for the derivation.

We reproduce Wei and Kirby (1995) with $d = 0.45$ m, $\delta = A/d = 0.1$, and a 450 m channel with (i) the present conservative staggered scheme, second order in space and time, and (ii) the HLLC scheme with fifth-order flux reconstruction and third-order RK time integration. In both cases, dispersive terms use second-order finite differences.
<figure><img src="image_20.png" /><figcaption>**Fig. 17.** Solitary wave test (Wei and Kirby, 1995): $\eta$ at $t = 150$ s with (a) the conservative staggered scheme; (b) with HLLC; both for $L/\Delta x = \{18, 36, 72, 144, 288\}$. (c) $L_\infty$ error vs. $L/\Delta x$ with fitted rates.</figcaption></figure>

Fig. 17 shows grid convergence to the analytical solitary-wave profile. HLLC exhibits a steeper pre-asymptotic rate ($p \approx 2.6$) and reaches a grid-converged profile by $L/\Delta x \approx 72$, whereas the conservative staggered scheme converges by $L/\Delta x \approx 144$ with $p \approx 1.7$. Beyond these thresholds, the apparent rates flatten because the reference solitary wave only approximates Nwogu's equations. The converged profiles agree well with the target, with a slightly overpredicted crest and weak trailing-edge dispersion, consistent with earlier studies (e.g., Kazolea and Delis (2013), Wei and Kirby (1995), Roeber et al. (2010)).

### A.3. Monochromatic wave propagation with $kd = \pi$

For completeness, we also present results for a shorter wave case ($kd = \pi$), which corresponds to the shortest resolvable wavelength supported by Nwogu's equations. The trends are consistent with those observed for $kd = 0.5$: the conservative staggered scheme exhibits a negative phase dispersion error, the HLLC scheme a positive phase

dispersion error, and the diffusive error remains negligible for the staggered scheme. The only differences appear in the magnitude of the phase shift and amplitude decay (see Fig. 18).

### A.4. Instability with the hybrid wave-breaking method

When using the hybrid wave-breaking approach with a grid size of $\Delta x = 0.05$ m, the simulation becomes unstable before reaching the full 2400 s runtime. The onset of the instability was observed at approximately 1245.1 s, where the transition from the dispersive to the hydrostatic solution at the breaking wave front ($x = 33$ m) induce spurious oscillations. Fig. 19 shows the evolution of the computed free surface just before the onset of the instability. As the shoaling wave begins to break at $t = 1244.9$ s, oscillations emerge at the interface between the breaking wave front, computed using a hydrostatic solution, and the downstream part of the wave. These oscillations continue to grow leading the breakdown of the numerical solution.
<figure><img src="image_21.png" /><figcaption>**Fig. 18.** Free-surface time series for monochromatic wave propagation with $kd = \pi$ and $a_0/d = 5 \times 10^{-3}$ at $x = 25\lambda_0$. A comparison between the present staggered scheme (second-order in space and time) and the HLLC scheme (third-order time integration and fifth-order reconstruction) for different grid sizes.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Free surface elevation at different time-stamps just before instability starts. The blue line represents the free surface computed using the full dispersive model, while the red lines highlight the regions identified as breaking grid cells, which are computed using the hyperbolic shallow water equation.</figcaption></figure>
# Data availability

The code that supports the findings of this study is openly available at: https://github.com/fzmihami/WavePropError.

## References

Abbott, M.B., Petersen, H.M., Skovgaard, O., 1978. On the numerical modelling of short waves in shallow water. J. Hydraul. Res. 16 (3), 173-204.

Barthelemy, X, Banner, ML, Peirson, WL, Fedele, F, Allis, M, Dias, F, 2018. On a unified breaking onset threshold for gravity waves in deep and intermediate depth water. J. Fluid Mech. 841, 463-488.

Beji, Serdar, Nadaoka, K., 1996. A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth. Ocean Eng. 23 (8), 691-704.

Bondehagen, Andreas, Roeber, Volker, Kalisch, Henrik, Buckley, Marc P, Streßer, Michael, Cysewski, Marius, Horstmann, Jochen, Bjørnestad, Maria, Ige, Olufemi E, Froysa, Hege G, et al., 2024. Wave-driven current and vortex patterns at an open beach: Insights from phase-resolving numerical computations and Lagrangian measurements. Coast. Eng. 193, 104591.

Boussinesq, Joseph, 1872. Théorie des ondes et des remous qui se propagent le long d'un canal rectangulaire horizontal, en commiquant au liquide contenu dans ce canal des vitesses sensiblement pareilles de la surface au fond. J. Mathématiques Pures et Appliquées 17, 55-108.

Choi, Young-Kwang, Shi, Fengyan, Malej, Matt, Smith, Jane M., 2018. Performance of various shock-capturing-type reconstruction schemes in the Boussinesq wave model, FUNWAVE-TVD. Ocean. Model. 131, 86-100.

Cienfuegos, Rodrigo, Barthélemy, Eric, Bonneton, Philippe, 2006. A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: model development and analysis. Internat. J. Numer. Methods Fluids 51 (11), 1217-1253.

Dingemans, M., 1973. Water Waves Over an Uneven Bottom; a Discussion of Long-Wave Equations. Technical report, Delft Hydraulics Laboratory.

Doyen, David, Gunawan, Putu Harry, 2014. An explicit staggered finite volume scheme for the Shallow Water Equations. In: Finite Volumes for Complex Applications VII-Methods and Theoretical Aspects. Springer, pp. 227-235.

Dutykh, Denys, Katsaounis, Theodoros, Mitsotakis, Dimitrios, 2011. Finite volume schemes for dispersive wave propagation and runup. J. Comput. Phys. 230 (8), 3035-3061.

Engquist, Björn, Osher, Stanley, 1981. One-sided difference approximations for nonlinear conservation laws. Math. Comp. 36 (154), 321-351.

Erduran, K.S., Ilic, S., Kutija, V., 2005. Hybrid finite-volume finite-difference scheme for the solution of Boussinesq equations. Internat. J. Numer. Methods Fluids 49 (11), 1213-1232.

Filippini, Andrea Gilberto, Bellec, Stevan, Colin, Mathieu, Ricchiuto, Mario, 2015. On the nonlinear behaviour of Boussinesq type models: Amplitude-velocity vs amplitude-flux forms. Coast. Eng. 99, 109-123.

Filippini, Andrea Gilberto, Kazolea, Maria, Ricchiuto, Mario, 2016. A flexible genuinely nonlinear approach for nonlinear wave propagation, breaking and run-up. J. Comput. Phys. 310, 381-417.

Fornberg, Bengt, 1998. A Practical Guide to Pseudospectral Methods, vol. 1, Cambridge University Press.

Fraccarollo, Luigi, Toro, Eleuterio F., 1995. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. J. Hydraul. Res. 33 (6), 843-864.

Gottlieb, Sigal, Shu, Chi-Wang, Tadmor, Eitan, 2001. Strong stability-preserving high-order time discretization methods. SIAM Rev. 43 (1), 89-112.

Gunawan, Harry Putu, 2015. Numerical simulation of Shallow Water Equations and related models (Ph.D. thesis). Paris Est.

Hamdi, S., Enright, W.H., Ouellet, Y., Schiesser, W.E., 2005. Method of lines solutions of the extended Boussinesq equations. J. Comput. Appl. Math. 183 (2), 327-342.

Harten, Amiram, Lax, Peter D., Leer, Bram van, 1983. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Rev. 25 (1), 35-61.

Jiang, Guang-Shan, Shu, Chi-Wang, 1996. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126 (1), 202-228.

Karambas, Theophanis V., Memos, Constantine D., 2009. Boussinesq model for weakly nonlinear fully dispersive water waves. J. Waterw. Port Coast. Ocean. Eng. 135 (5), 187-199.

Kazolea, M., Delis, A.I., 2013. A well-balanced shock-capturing hybrid finite volume-finite difference numerical scheme for extended 1D Boussinesq models. Appl. Numer. Math. 67, 167-186.

Kazolea, Maria, Ricchiuto, Mario, 2018. On wave breaking for Boussinesq-type models. Ocean. Model. 123, 16-39.

Kennedy, Andrew B, Chen, Qin, Kirby, James T, Dalrymple, Robert A, 2000. Boussinesq modeling of wave transformation, breaking, and runup. I: 1D. J. Waterw. Port Coast. Ocean. Eng. 126 (1), 39-47.

Kim, Kyu Hong, Kim, Chongam, 2005. Accurate, efficient and monotonic numerical methods for multi-dimensional compressible flows: Part II: Multi-dimensional limiting process. J. Comput. Phys. 208 (2), 570-615.

Kim, Dae-Hong, Lynett, Patrick J., Socolofsky, Scott A., 2009. A depth-integrated model for weakly dispersive, turbulent, and rotational fluid flows. Ocean. Model. 27 (3-4), 198-214.

Kirby, James T., 2016. Boussinesq models and their application to coastal processes across a wide range of scales. J. Waterw. Port Coast. Ocean. Eng. 142 (6), 03116005.

Kirby, James T., 2020. A new instability for Boussinesq-type equations. J. Fluid Mech. 894, F1.

Kolmogorov, Andrej Nikolaevich, 1941. Equations of turbulent motion in an incompressible fluid. In: Dokl. Akad. Nauk SSSR, vol. 30, pp. 299-303.

Kurganov, Alexander, Petrova, Guergana, 2007. A second-order well-balanced positivity preserving central-upwind scheme for the Saint-Venant system. Commun. Math. Sci. 5 (1), 133-160.

Larsen, Jesper, Dancy, Henry, 1983. Open boundaries in short wave simulations—a new approach. Coast. Eng. 7 (3), 285-297.

Le Méhauté, Bernard, 1976. An Introduction to Hydrodynamics and Water Waves. Springer Science & Business Media.

Lynett, P., Liu, P.L.F., Sitanggang, K.I., Kim, D.H., 2002. Modeling Wave Generation, Evolution, and Interaction with Depthintegrated, Dispersive Wave Equations COULWAVE Code Manual. Technical report, Cornell University.

Lynett, Patrick J, Swigler, David, Son, Sangyoung, Bryant, Duncan, Socolofsky, Scott, 2010. Experimental study of solitary wave evolution over a 3D shallow shelf. Coast. Eng. Proc. x (32), 1-1.

Madsen, P.A., Schäffer, H.A., 1998. Higher-order Boussinesq-type equations for surface gravity waves: derivation and analysis. Philos. Trans. R. Soc. Lond. Ser. A. Math. Phys. Eng. Sci. 356 (1749), 3123-3181.

Madsen, Per A., Simonsen, Hans Jacob, Pan, Cun-Hong, 2005. Numerical simulation of tidal bores and hydraulic jumps. Coast. Eng. 52 (5), 409-433.

Madsen, Per A., Sørensen, Ole R., 1992. A new form of the Boussinesq equations with improved linear dispersion characteristics. Part 2. A slowly-varying bathymetry. Coast. Eng. 18 (3-4), 183-204.

Madsen, Per A., Sørensen, O.R., Schäffer, H.A., 1997. Surf zone dynamics simulated by a Boussinesq type model. Part I. Model description and cross-shore motion of regular waves. Coast. Eng. 32 (4), 255-287.

Mihami, Fatima Zahra, 2025. WavePropError: A Hybrid Python-C++ Solver for the One-Dimensional Nwogu Equations. https://github.com/fzmihami/WavePropError. GitHub repository.

Mihami, Fatima-Zahra, Roeber, Volker, Morichon, Denis, 2022. Efficient numerical computations of long-wave run-up and their sensitivity to grid nesting. Water Waves 4 (3), 517-548.

Nwogu, Okey, 1993. Alternative form of Boussinesq equations for nearshore wave propagation. J. Waterw. Port, Coast. Ocean. Eng. 119 (6), 618-638.

Nwogu, Okey, George, 1996. Numerical prediction of breaking waves and currents with a Boussinesq model. In: Proceedings of the 25th International Conference on Coastal Engineering. ICCE, http://dx.doi.org/10.9753/icce.v25.

Peregrine, D. Howell, 1967. Long waves on a beach. J. Fluid Mech. 27 (4), 815-827.

Pope, S.B., 2000. Turbulent Flows. Journal of Turbulence, vol. 1, Cambridge University Press, Cambridge, p. 771. http://dx.doi.org/10.1088/1468-5248/1/1/702.

Prandtl, Ludwig, et al., 1945. Über ein neues formelsystem für die ausgebildete turbulenz. Nacr. Akad. Wiss. Gottingen, Math-Phys. Kl 6-19.

Roe, Philip L., 1997. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 135 (2), 250-258.

Roeber, Volker, Cheung, Kwok Fai, 2012. Boussinesq-type model for energetic breaking waves in fringing reef environments. Coast. Eng. 70, 1-20.

Roeber, Volker, Cheung, Kwok Fai, Kobayashi, Marcelo H., 2010. Shock-capturing Boussinesq-type model for nearshore wave processes. Coast. Eng. 57 (4), 407-423.

Roelvink, Dano, McCall, Robert, Mehvar, Seyedabdolhossein, Nederhoff, Kees, Dastgheib, Ali, 2018. Improving predictions of swash dynamics in XBeach: The role of groupiness and incident-band runup. Coast. Eng. 134, 103-123.

Shi, Fengyan, Kirby, James T, Harris, Jeffrey C, Geiman, Joseph D, Grilli, Stephan T, 2012. A high-order adaptive time-stepping TVD solver for Boussinesq modeling of breaking waves and coastal inundation. Ocean. Model. 43, 36-51.

Shi, Fengyan, Kirby, James T, Tehranirad, Babak, Harris, Jeffrey C, Grilli, ST, 2011. FUNWAVE-TVD Fully Nonlinear Boussinesq Wave Model with TVD Solver Documentation and User's Manual. Cent. Appl. Coast. Res. Univ. Del. Newark USA Res. Rep. No CACR-11-04 2, 1.

Sørensen, Ole R, Schäffer, Hemming A, Sørensen, Lars S, 2004. Boussinesq-type modelling using an unstructured finite element technique. Coast. Eng. 50 (4), 181-198.

Stelling, Guss S., Duinmeijer, S.P. Alex, 2003. A staggered conservative scheme for every Froude number in rapidly varied shallow water flows. Internat. J. Numer. Methods Fluids 43 (12), 1329-1354.

Stelling, Guus, Zijlema, Marcel, 2003. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation. Internat. J. Numer. Methods Fluids 43 (1), 1-23.

Suzuki, Tomohiro, Altomare, Corrado, Veale, William, Verwaest, Toon, Trouw, Koen, Troch, Peter, Zijlema, Marcel, 2017. Efficient and robust wave overtopping estimation for impermeable coastal structures in shallow foreshores using SWASH. Coast. Eng. 122, 108-123.
Swigler, David Townley, 2010. Laboratory Study Investigating the Three-Dimensional Turbulence and Kinematic Properties Associated with a Breaking Solitary Wave (Ph.D. thesis). Texas A & M University.

Tissier, M., Bonneton, P., Marche, F., Chazel, F., Lannes, D., 2012. A new approach to handle wave breaking in fully non-linear Boussinesq models. Coast. Eng. 67, 54–66.

Tonelli, Mara, Petti, Marco, 2009. Hybrid finite volume–finite difference scheme for 2DH improved Boussinesq equations. Coast. Eng. 56 (5–6), 609–620.

Tonelli, Mara, Petti, Marco, 2011. Simulation of wave breaking over complex bathymetries by a Boussinesq model. J. Hydraul. Res. 49 (4), 473–486.

Toro, Eleuterio, 2001. Shock-capturing methods for free-surface shallow flows, vol. 868, Wiley New York.

Toro, E.F., Chakraborty, A., 1994. The development of a Riemann solver for the steady supersonic Euler equations. Aeronaut. J. 98 (979), 325–339.

Torres, Marissa J, Lam, Michael-Angelo Y-H, Malej, Matt, et al., 2022. Practical Guidance for Numerical Modeling in FUNWAVE-TVD. Technical report, Engineer Research and Development Center (US).

Varing, Audrey, Filipot, Jean-Francois, Grilli, Stephan, Duarte, Rui, Roeber, Volker, Yates, Marissa, 2021. A new definition of the kinematic breaking onset criterion validated with solitary and quasi-regular waves in shallow water. Coast. Eng. 164, 103755.

Wei, Ge, Kirby, James T., 1995. Time-dependent numerical code for extended Boussinesq equations. J. Waterw. Port Coast. Ocean. Eng. 121 (5), 251–261.

Wei, Ge, Kirby, James T., Sinha, Amar, 1999. Generation of waves in Boussinesq models using a source function method. Coast. Eng. 36 (4), 271–299.

Yamamoto, S, Daiguji, H, 1993. Higher-order-accurate upwind schemes for solving the compressible Euler and Navier-Stokes equations. Comput. & Fluids 22 (2–3), 259–270.

Yamazaki, Yoshiki, Cheung, Kwok Fai, Kowalik, Zygmunt, 2011. Depth-integrated, non-hydrostatic model with grid nesting for tsunami generation, propagation, and run-up. Internat. J. Numer. Methods Fluids 67 (12), 2081–2107.

Yamazaki, Yoshiki, Cheung, Kwok Fai, Kowalik, Z. Lay, T, Pawlak, G, 2012. Neowave. In: Proceedings and Results of the 2011 NTHMP Model Benchmarking Workshop, Boulder: US Department of Commerce/NOAA/NTHMP (NOAA Special Report). pp. 239–302.

Yamazaki, Yoshiki, Kowalik, Zygmunt, Cheung, Kwok Fai, 2009. Depth-integrated, non-hydrostatic model for wave breaking and run-up. Internat. J. Numer. Methods Fluids 61 (5), 473–497.

Yuan, Ye, Shi, Fengyan, Kirby, James T., Yu, Fujiang, 2020. FUNWAVE-GPU: Multiple-GPU Acceleration of a Boussinesq-Type Wave Model. J. Adv. Model. Earth Syst. 12 (5), e2019MS001957.

Zelt, J.A., 1991. The run-up of nonbreaking and breaking solitary waves. Coast. Eng. 15, 205–246.

Zhang, Mingliang, Wu, W.M., 2011. A two dimensional hydrodynamic and sediment transport model for dam break based on finite volume method with quadtree grid. Appl. Ocean Res. 33 (4), 297–308.

Zhao, Kuifeng, Wang, Yufei, Liu, Philip L.-F., 2024. A guide for selecting periodic water wave theories-Le Méhauté (1976)'s graph revisited. Coast. Eng. 188, 104432.

Zhou, J.G., Stansby, P.K., 1999. 2D shallow water flow model for the hydraulic jump. Internat. J. Numer. Methods Fluids 29 (4), 375–387.

Zhou, Quan, Zhan, Jiemin, Li, Yoksheung, 2016. High-order finite volume WENO schemes for Boussinesq modelling of nearshore wave processes. J. Hydraul. Res. 54 (6), 646–662.

Zijlema, Marcel, 2019. The role of the Rankine-Hugoniot relations in staggered finite difference schemes for the Shallow Water Equations. Comput. & Fluids 192, 104274.

Zijlema, Marcel, Stelling, Guus, Smit, Pieter, 2011. SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coast. Eng. 58 (10), 992–1012.