2. Generalities on the depth-averaged fields

In this section, we introduce the depth-averaging procedure, and apply it to the local velocity and vorticity fields. This allows us to draw some general preliminary observations that will be useful later on.

For the sake of simplicity, we adopt the following notation for the two- and three-dimensional fields. Let us denote by ${\boldsymbol {u}} = ( u, v, w )$ the three-dimensional velocity field. Then we indicate by ${\overline {\boldsymbol {u}}} = ( u, v )$ the two-dimensional velocity field on the $(x,y)$ -plane. Generally, the same notation is adopted to indicate the projection on the $(x,y)$ -plane of a generic three-dimensional field. Similarly, we denote by $\overline {\nabla }$ the gradient operator on the $(x,y)$ -plane, while $\nabla$ is the gradient in three dimensions. Finally, the quantities evaluated at the free surface and along the seabed are denoted by the subscripts $F$ and $B$ . In particular, the elevation of the free surface is indicated by $\eta (x,y,t)$ , while the seabed position is given by $z=-h(x,y,t)$ . Accordingly, the total water depth is $d = \eta + h$ . The origin of the Cartesian frame of reference is placed at the position of the undisturbed free surface, with the vertical axis pointing upwards.

The components of the vorticity field $\boldsymbol {\omega } = \nabla \times {\boldsymbol {u}}$ are

(2.1) \begin{eqnarray} \omega _1 = \frac {\partial w}{\partial y} - \frac {\partial v}{\partial z} , \quad \omega _2 = \frac {\partial u}{\partial z} - \frac {\partial w}{\partial x} , \quad \omega _3 = \frac {\partial v}{\partial x} -\frac {\partial u}{\partial y} . \end{eqnarray}

Integrating the first two relations from $z$ to the free surface $\eta$ along the vertical direction, using the Leibniz rule and rearranging, we find

(2.2) \begin{align} u &= - \frac {\partial \Upsilon }{\partial x} + u_F + w_F \, \frac {\partial \eta }{\partial x} - {\int _{z}^{\eta }} \omega _2 \, {\textrm {d}}\zeta , \end{align}

(2.3) \begin{align} v &= - \frac {\partial \Upsilon }{\partial y} + v_F + w_F \, \frac {\partial \eta }{\partial y} + {\int _{z}^{\eta }} \omega _1 \, {\textrm {d}}\zeta , \end{align}

where the variable $\Upsilon$ is defined by the expression

(2.4) \begin{eqnarray} \Upsilon = {\int _{z}^{\eta }} w \,{\textrm {d}}\zeta \quad \Rightarrow \quad w = - \frac {\partial \Upsilon }{\partial z} . \end{eqnarray}

We observe that the function $\Upsilon$ behaves as a velocity potential for $\boldsymbol {u}$ , whereas the remaining contributions, at this point, do not yet have a clear interpretation. In order to provide a physical characterization, we decompose them in a depth-averaged part and a rotational deviation. This leads to the definition of the generalized mass flux ${\boldsymbol {M}} = (M_1 , M_2 )$ as

(2.5) \begin{align} \frac {M_1}{d} &= u_F + w_F \, \frac {\partial \eta }{\partial x} - \frac {1}{d} {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \omega _2 \, {\textrm {d}}\zeta , \end{align}

(2.6) \begin{align} \frac {M_2}{d} &= v_F + w_F \, \frac {\partial \eta }{\partial y} + \frac {1}{d} {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \omega _1 \, {\textrm {d}}\zeta , \end{align}

so that (2.2) and (2.3) can be rewritten as

(2.7) \begin{eqnarray} {\overline {\boldsymbol {u}}} = - \overline {\nabla } \Upsilon + \frac {{\boldsymbol {M}}}{d} + {\boldsymbol {R}} , \end{eqnarray}

where the rotational deviation is given by

(2.8) \begin{eqnarray} {\boldsymbol {R}} = \left (- {\int _{z}^{\eta }} \omega _2 \, {\textrm {d}}\zeta + \frac {1}{d} {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \omega _2 \, {\textrm {d}}\zeta {\int _{z}^{\eta }} \omega _1 \, {\textrm {d}}\zeta - \frac {1}{d} {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \omega _1 \, {\textrm {d}}\zeta\right ) . \end{eqnarray}

As shown in § 3.2, the field $\boldsymbol {M}$ naturally appears when depth averaging is applied to the Navier–Stokes equations, and as pointed out in Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017), it includes those linear contributions that are associated with the phenomenon of wave dispersion. This motivates the wording ‘generalized’ adopted for it. Accordingly, we refer to the field ${\boldsymbol {M}}/d$ as the generalized velocity field. Finally, integrating (2.7) over the fluid depth, we obtain

(2.9) \begin{eqnarray} \textrm{d} {\boldsymbol {U}} = - {\int _{-h}^{\eta }} \overline {\nabla } \Upsilon \, {\textrm {d}z} + {{\boldsymbol {M}}} , \end{eqnarray}

where $\boldsymbol {U} = (U,V)$ is the depth-averaged field on the $(x,y)$ -plane, i.e.

(2.10) \begin{eqnarray} \boldsymbol {U} = \frac {1}{d} {\int _{-h}^{\eta }} {\overline {\boldsymbol {u}}} \, {\textrm {d}z} . \end{eqnarray}

Accordingly, we denote by ${\delta \boldsymbol {u}} = {\overline {\boldsymbol {u}}} - \boldsymbol {U}$ the deviation of the two-dimensional velocity field from the depth-averaged quantities.

Before proceeding, it is worth noting that some kinds of ‘generalized velocities’ also appear in many models for coastal dynamics, even though there is not always a full awareness of it. For example, in Wei et al. (Reference Wei, Kirby, Grilli and Subramanya1995), Erduran, Ilic & Kutija (Reference Erduran, Ilic and Kutija2005) and Tonelli & Petti (Reference Tonelli and Petti2012), the variable in the momentum equation that is integrated in time is not the depth-averaged velocity, rather an expression containing both $\boldsymbol {U}$ and its spatial derivatives (the latter extracted from the fluxes of the momentum equation). Such a variable can be regarded as a generalized velocity field due to the similarities with the derivation of $\boldsymbol {M}/d$ in § 3.2. In other models (Nwogu Reference Nwogu1993; Chen et al. Reference Chen, Kirby, Dalrymple, Kennedy and Chawla2000; Kennedy et al. Reference Kennedy, Chen, Kirby and Dalrymple2000), in addition to $\boldsymbol {U}$ , a further velocity term is defined by evaluating $\overline {\boldsymbol {u}}$ at a certain reference quote over the fluid depth. This strategy is applied to improve the description of the dispersive effects in wave propagation. Again, this procedure shows close analogies with the physical interpretation of $\boldsymbol {M}/d$ as a term consisting of velocity components along the free surface; see (2.5)–(2.6).

As a final step, we highlight some aspects related to the decomposition provided in this section. By a slight abuse of notation, we introduce the curl operator in two dimensions as $\overline {\nabla } \times \boldsymbol {f} = \partial f_2/\partial x - \partial f_1/\partial y$ , where $\boldsymbol {f} = (f_1, f_2)$ is a generic vector. Hence applying it to (2.7), we find

(2.11) \begin{eqnarray} \omega _3 = \overline {\nabla } \times \left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R} \right ) . \end{eqnarray}

The above result, already obtained in Antuono et al. (Reference Antuono, Lucarelli, Bardazzi and Brocchini2022), shows that $\omega _3$ can be derived from the remaining components of the vorticity field. This is a consequence of the fact that $\boldsymbol {\omega }$ is a null-divergence field. Furthermore, (2.11) suggests a close relation between the field $\boldsymbol {M}/d$ and the vorticity in the fluid. This aspect will be addressed in detail in the following sections.

2.1. Rotational and gradient fields

Let us introduce the normal unit vectors to the free surface and the seabed. Since we require them to point out of the fluid domain, they read

(2.12) \begin{eqnarray} \boldsymbol {n}_F = \left ( -\frac {\partial \eta }{\partial x}, -\frac {\partial \eta }{\partial y}, 1 \right ) / N_F, \quad \boldsymbol {n}_B = \left ( -\frac {\partial h}{\partial x}, -\frac {\partial h}{\partial y}, -1 \right ) / N_B, \end{eqnarray}

where $N_F^2 = 1 + \| \overline {\nabla } \eta \|^2$ and $N_B^2 = 1 + \| \overline {\nabla } h \|^2$ . We also define $\tilde {\boldsymbol {n}}_F = N_F \boldsymbol {n}_F$ and $\boldsymbol {\tilde {n}}_B = N_B \,\boldsymbol {n}_B$ , since these vectors will be used extensively for the derivation of useful expressions. Note that $\boldsymbol {\tilde {n}}_F$ and $\boldsymbol {\tilde {n}}_B$ are not unit vectors.

In particular, as shown in Appendix A, the following relation holds true:

(2.13) \begin{eqnarray} \overline {\nabla } \times (\boldsymbol {R}_F - \boldsymbol {R}_B) = \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F + \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B , \end{eqnarray}

where

(2.14) \begin{eqnarray} \boldsymbol {R}_F - \boldsymbol {R}_B = \left ({\int _{-h}^{\eta }} \omega _2 \, {\textrm {d}z} , - {\int _{-h}^{\eta }} \omega _1 \, {\textrm {d}z} \right ) . \end{eqnarray}

In more detail, the following expressions between the vorticity along the interfaces and the generalized velocity $\boldsymbol {M}/d$ are derived:

(2.15) \begin{align} \overline {\nabla } \times \left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right ) &= { \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F }, \end{align}

(2.16) \begin{align} \overline {\nabla } \times \left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right ) &= { - \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B }. \end{align}

The above relations are trivially satisfied when unidirectional wave propagation occurs. Conversely, for a generic wave motion, they suggest a special structure for the equations that describe the evolution of the generalized velocity. For example, the term $\boldsymbol {M}/d$ is expected to be a gradient field when the flow is irrotational. The same reasoning stands for the fields $( \boldsymbol {M}/d + \boldsymbol {R}_F )$ and $( \boldsymbol {M}/d + \boldsymbol {R}_B )$ if the boundary terms $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ are null. Apart from the above considerations, (2.15) and (2.16) suggest that these boundary terms may play an important role in the characterization of the depth-averaged models. These relevant aspects will be addressed in the following sections where the equations of motions are studied in the framework of the depth-averaging approach. Incidentally, we observe that, differently from $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ , $\boldsymbol {R}_F$ and $\boldsymbol {R}_B$ do not represent boundary conditions, since they are defined through integrals over the water depth and therefore contain vortical contributions coming from both the interfaces (namely, free surface and seabed).

Finally, we highlight that the expressions and results derived up to here represent general properties of the depth-averaged field, and do not depend on any specific dynamic law.

3. Governing equations and depth averaging

In this section, we apply the depth-averaging procedure to the equations for incompressible Newtonian fluids. The analysis will be developed by keeping in mind the results of the previous section and highlighting the presence of terms depending on $( \boldsymbol {M}/d + \boldsymbol {R}_F )$ , $( \boldsymbol {M}/d + \boldsymbol {R}_B )$ and the boundary expressions $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ .

Incidentally, we clarify that in depth-averaged models, the interfaces (i.e. the free surface and the seabed) are assumed to be single-valued in order to avoid ambiguities during the integration along the vertical direction. This implies that the vortical contributions generated by wave breaking (which generally involve free-surface overturning and/or fragmentation) have to be modelled separately through proper closures. This approach also stands for the vorticity generated at the seabed, since the spatial resolution used in the numerical implementation is generally too coarse to describe the boundary layer. In practice, all the dynamics of the rotational component of the velocity field has to be modelled, apart from the main depth-averaged equations. In the present work, we do not deal with the latter part, which represents a distinct line of research, and assume that a proper model is available for the vorticity field. Similarly, the same assumption is adopted for the turbulent terms whose closure is not discussed here.

To simplify the treatise as much as possible, we write the Navier–Stokes equations in the following compact form:

(3.1) \begin{eqnarray} \left \{\! \begin{array}{l} \nabla \cdot {\boldsymbol {u}} = 0,\\ \displaystyle \frac {\partial {\boldsymbol {u}}}{\partial t} = - \nabla P + \boldsymbol {g} - \boldsymbol {q} + \boldsymbol {\theta }, \end{array} \right . \end{eqnarray}

where $P = p + \| {\boldsymbol {u}} \|^2/2$ , $p$ is the pressure field (divided by the density), $\boldsymbol {g}$ is the gravity acceleration, $\boldsymbol {q} = \boldsymbol {\omega } \times {\boldsymbol {u}}$ is the Lamb vector, and $\boldsymbol {\theta } = \nabla \cdot \mathbb {V}$ , where $\mathbb {V}$ is the viscous stress tensor. The latter terms also includes possible turbulent contributions. The kinematic boundary conditions are

(3.2) \begin{eqnarray} w_F = \frac {\partial \eta }{\partial t} + {\overline {\boldsymbol {u}}}_F \cdot \overline {\nabla } \eta , \quad w_B = - \frac {\partial h}{\partial t} - {\overline {\boldsymbol {u}}}_B \cdot \overline {\nabla } h , \end{eqnarray}

while the dynamics boundary condition at the free surface is

(3.3) \begin{eqnarray} p_F = \left ( \boldsymbol {n} \cdot \mathbb {V} \cdot \boldsymbol {n} \right )_F . \end{eqnarray}

Here it is also possible to include the action of surface tension. Integrating the continuity equation over the water depth and using the kinematic conditions, we find:

(3.4) \begin{eqnarray} \frac {\partial d}{\partial t} \, + \, \overline {\nabla } \cdot \left ( \boldsymbol {U}\, d\right ) \, = \, 0 \,. \end{eqnarray}

From a certain point of view, (3.4) represents the depth-averaged counterpart of the kinematic conditions along the free surface and the seabed. Substituting (2.2)–(2.4) in the continuity equation, we finally find the governing equation for $\Upsilon$ ,

(3.5) \begin{eqnarray} \nabla \Upsilon = \nabla \cdot \left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R} \right ) , \end{eqnarray}

along with the following boundary conditions at the free surface and at the seabed:

(3.6) \begin{eqnarray} \Upsilon _F = 0 , \quad \frac {\partial \Upsilon }{\partial n} \bigg |_B = - \left [\frac {\partial h}{\partial t} + \left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right ) \cdot \overline {\nabla } h \right ]/N_B . \end{eqnarray}

The Neumann condition is equivalent to the kinematic condition along the seabed, and its derivation is detailed in Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017). Equation (3.5) for $\Upsilon$ is solved in place of the Poisson equation for the pressure field that is obtained by taking the divergence of the momentum equation in (3.1). In particular, the solution of the Poisson equation for $\Upsilon$ is used to recover the dynamics along the vertical direction relying on the knowledge of the fields $\boldsymbol {M}/d$ and $\boldsymbol {R}$ . As illustrated in the following sections, the governing equation for $\boldsymbol {M}/d$ is obtained from the momentum equations through the application of the depth-averaging procedure.

3.1. The equation for the dynamic pressure field

Integrating the third component of the momentum equation over the vertical interval $[z,\eta ]$ , we obtain

(3.7) \begin{eqnarray} {\int _{z}^{\eta }} \frac {\partial w}{\partial t} \, {\textrm {d}}\zeta = - P_F + P - g ( \eta -z ) - {\int _{z}^{\eta }} q_3 \, {\textrm {d}}\zeta + {\int _{z}^{\eta }} \theta _3 \, {\textrm {d}}\zeta, \end{eqnarray}

from which we find

(3.8) \begin{eqnarray} { P = \frac {\partial \Upsilon }{\partial t} + P_F - w_F \, \frac {\partial \eta }{\partial t} + g ( \eta -z ) + {\int _{z}^{\eta }} q_3 \, {\textrm {d}}\zeta - {\int _{z}^{\eta }} \theta _3 \, {\textrm {d}}\zeta .} \end{eqnarray}

Similarly, integrating over $[-h, \eta ]$ , we obtain

(3.9) \begin{eqnarray} { P_B = \frac {\partial \Upsilon _B }{\partial t} + P_F - w_F \, \frac {\partial \eta }{\partial t} - w_B \, \frac {\partial h}{\partial t} + g d + {\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} - {\int _{-h}^{\eta }} \theta _3 \, {\textrm {d}z} .} \end{eqnarray}

Finally, we recall that

(3.10) \begin{eqnarray} P_F = p_F + \frac {\| {\boldsymbol {u}}_F \|^2}{2} , \end{eqnarray}

where $p_F$ is given by (3.3). The above expressions are used to eliminate the explicit dependence on the pressure field from the remaining components of the momentum equation.

Before proceeding to the analysis, we highlight that the derivation of many Boussinesq and non-hydrostatic models can be basically traced back to approximations applied to the dynamics pressure field (e.g. Nwogu Reference Nwogu1993; Kim et al. Reference Kim, Lynett and Socolofsky2009; Donahue et al. Reference Donahue, Zhang, Kennedy, Westerink, Panda and Dawson2015; Kazakova & Richard Reference Kazakova and Richard2019) and/or to the Laplace equation for the velocity potential $\phi$ (e.g. Wei et al. Reference Wei, Kirby, Grilli and Subramanya1995; Gobbi et al. Reference Gobbi, Kirby and Wei2000; Raoult et al. Reference Raoult, Benoit and Yates2016). In the current framework, the latter case is equivalent to approximating the solution of the Poisson equation for $\Upsilon$ , namely (3.5), whereas the former approach corresponds to approximations of (3.8)–(3.10). Regardless of the accuracy of the above strategies, an insidious issue lies in the possible loss of general physical properties associated with the wave dynamics. For this reason, in the following analysis we do not consider any of the above approximations, and derive a general structure that stands as a reference framework for both Boussinesq and non-hydrostatic models.

3.2. The equation for the generalized velocity field

Let us consider the momentum in the $x$ -direction. Integrating over the water depth and using the Leibniz rule, we find

(3.11) \begin{align} \frac {\partial ( {U} d )}{\partial t} - u_F \, \frac {\partial \eta }{\partial t} - u_B \, \frac {\partial h}{\partial t} &= - \frac {\partial }{\partial x} \left ( {\int _{-h}^{\eta }} P \, {\textrm {d}z} \right ) + P_F \, \frac {\partial \eta }{\partial x} + P_B \, \frac {\partial h}{\partial x} \nonumber \\ &\quad {}- {\int _{-h}^{\eta }} q_1 \, {\textrm {d}z} + {\int _{-h}^{\eta }} \theta _1 \, {\textrm {d}z} . \end{align}

We first focus on the integral of the pressure field. Substituting (3.8) and using the condition $\Upsilon _F=0$ and the Leibniz rule, we find

(3.12) \begin{eqnarray} {\int _{-h}^{\eta }} P \, {\textrm {d}z} &=& \frac {\partial }{\partial t}\left ( {\int _{-h}^{\eta }} \Upsilon \, {\textrm {d}z} \right ) - \Upsilon _B \, \frac {\partial h}{\partial t} + d \left ( P_F - w_F\, \frac {\partial \eta }{\partial t} \right ) + g \, \frac {d^2}{2} \\ && {}+{\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} q_3 \, {\textrm {d}}\zeta - {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \theta _3 \, {\textrm {d}}\zeta .\nonumber \end{eqnarray}

Substituting (3.12) and (3.9) in (3.11) and rearranging, we obtain

(3.13) \begin{eqnarray} \frac {\partial }{\partial t}\left [{U} d + \frac {\partial }{\partial x}\left ( {\int _{-h}^{\eta }} \Upsilon \,{\textrm {d}z} \right ) - \Upsilon _B \, \frac {\partial h}{\partial x} \right ] = u_F \, \frac {\partial \eta }{\partial t} + u_B \, \frac {\partial h}{\partial t} + \frac {\partial \Upsilon _B}{\partial x} \, \frac {\partial h}{\partial t} - g d \, \frac {\partial \eta }{\partial x} \nonumber \\ - \frac {\partial }{\partial x}\left [d \left (P_F - w_F \, \frac {\partial \eta }{\partial t} \right ) \right ] + P_F \, \frac {\partial d}{\partial x} - \left ( w_F \, \frac {\partial \eta }{\partial t} + w_B \, \frac {\partial h}{\partial t} \right ) \frac {\partial h}{\partial x} + G_1\; , \end{eqnarray}

where

(3.14) \begin{eqnarray} G_1 &=& - \frac {\partial }{\partial x}\left ({\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} q_3 \, {\textrm {d}}\zeta - {\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \theta _3 \, {\textrm {d}}\zeta \right ) \nonumber \\ &&{}+ \left ( {\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} - {\int _{-h}^{\eta }} \theta _3 \, {\textrm {d}z} \right ) \frac {\partial h}{\partial x} - {\int _{-h}^{\eta }} q_1 \, {\textrm {d}z} + {\int _{-h}^{\eta }} \theta _1 \, {\textrm {d}z} . \quad \end{eqnarray}

Hereinafter, we introduce the two-dimensional vector $\boldsymbol {G}= (G_1, G_2)$ , where $G_2$ is the component along the $y$ -direction analogous to $G_1$ . Note that (3.13) represents a cornerstone in the derivation of the models, since it provides the passage of certain time derivatives out of the momentum flux, and motivates the definition of the generalized mass flux $\boldsymbol {M}$ . Indeed, using (2.9) and rearranging, we find

(3.15) \begin{eqnarray} \frac {\partial M_1 }{\partial t} &=& d \, \frac {\partial }{\partial x}\left (w_F \, \frac {\partial \eta }{\partial t} - P_F - g \eta \right ) + \left (u_F + w_F \, \frac {\partial \eta }{\partial x} \right ) \frac {\partial \eta }{\partial t} \nonumber \\ && {}+\left (u_B - w_B \, \frac {\partial h}{\partial x} + \frac {\partial \Upsilon _B }{\partial x} \right ) \frac {\partial h}{\partial t} + G_1 . \end{eqnarray}

Now, let us focus on (2.7). Using (2.4) for $\Upsilon$ , we simplify some expressions as

(3.16) \begin{eqnarray} u_F + w_F \, \frac {\partial \eta }{\partial x} = \frac {{M}_1}{d} + \left ( {R}_F \right )_1 , \quad u_B - w_B \, \frac {\partial h}{\partial x} + \frac {\partial \Upsilon _B}{\partial x} = \frac {{M}_1}{d} + \left ( {R}_B \right )_1 , \end{eqnarray}

where the condition $\Upsilon _F = 0$ is applied in the first relation. Substituting these relations in the momentum equation, dividing by the total water depth and rearranging, we obtain a compact form that in vectorial notation reads as

(3.17) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} \right ) = \overline {\nabla } \left (w_F \, \frac {\partial \eta }{\partial t} - P_F - g \eta \right ) + \frac {1}{d} \left ( \boldsymbol {R}_F \, \frac {\partial \eta }{\partial t} + \boldsymbol {R}_B \, \frac {\partial h}{\partial t} + \boldsymbol {G} \right ) . \end{eqnarray}

The above expression is compatible with the condition $\overline {\nabla } \times (\boldsymbol {M}/d) = 0$ described in Antuono et al. (Reference Antuono, Valenza, Lugni and Colicchio2019) for inviscid and irrotational laminar flows (i.e. when $\boldsymbol {R}=0$ and $\boldsymbol {G}=0$ ). Besides this, it is possible to provide a deeper inspection of the above structure. We observe that

(3.18) \begin{eqnarray} \boldsymbol {R}_F \, \frac {\partial \eta }{\partial t} + \boldsymbol {R}_B \, \frac {\partial h}{\partial t} &=& \boldsymbol {R}_F \, \frac {\partial d}{\partial t} - ( \boldsymbol {R}_F - \boldsymbol {R}_B ) \, \frac {\partial h}{\partial t} \nonumber \\ &=& \frac {\partial (\boldsymbol {R}_F d)}{\partial t} - d \, \frac {\partial \boldsymbol {R}_F }{\partial t} - ( \boldsymbol {R}_F - \boldsymbol {R}_B ) \, \frac {\partial h}{\partial t}, \end{eqnarray}

and (3.17) becomes

(3.19) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right ) \, &=& \, \overline {\nabla } \left (\, w_F \, \frac {\partial \eta }{\partial t} - P_F - g\eta \,\right ) \nonumber \\ &&{} + \frac {1}{d} \left [\, \frac {\partial (\boldsymbol {R}_F \, d)}{\partial t} - ( \boldsymbol {R}_F \!-\! \boldsymbol {R}_B ) \, \frac {\partial h}{\partial t} + \boldsymbol {G} \,\right ] \,. \quad \end{eqnarray}

In the next subsection, we show that the last term in the right-hand side can be recast in the form of vortical and viscous boundary contributions at the free surface. This result represents a keystone in the present reformulation of the depth-averaged equations, since it allows for a drastic simplification of the equations themselves, and for the identification of the intrinsic boundary conditions along this interface.

3.3. The equation for the vorticity field

The equation for the vorticity is obtained by applying the curl operator to the momentum equation in the second expression of the system (3.1). We obtain

(3.20) \begin{eqnarray} \frac {\partial \boldsymbol {\omega }}{\partial t} = - \nabla \times \boldsymbol {q} + \nabla \times \boldsymbol {\theta } . \end{eqnarray}

As shown in § B.1, this equation along with the kinematic condition at the free surface allows one to obtain the result

(3.21) \begin{eqnarray} \frac {1}{d} \left [\frac {\partial (\boldsymbol {R}_F d)}{\partial t} - ( \boldsymbol {R}_F - \boldsymbol {R}_B ) \, \frac {\partial h}{\partial t} + \boldsymbol {G} \right ] = {\overline {\boldsymbol {u}}}^{\perp }_F \left ( \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F \right ) + \left [\overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F \,\overline {\nabla } \eta \right ], \end{eqnarray}

where ${\overline {\boldsymbol {u}}}^\perp = (v, - u)$ . The above expression allows us to simplify (3.19) as

(3.22) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right ) = \overline {\nabla } \left ( w_F \, \frac {\partial \eta }{\partial t} - P_F - g \eta \right ) + {\overline {\boldsymbol {u}}}^{\perp }_F \left ( \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F \right ) + \left [\overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F\, \overline {\nabla } \eta \right ] .\qquad \end{eqnarray}

This formulation is important since it brings to light the significant boundary conditions along the free surface in an explicit form. We observe that in addition to the term proportional $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ , there is a further boundary contribution that originates from the viscous term $\boldsymbol {\theta }$ .

One interesting aspect is that it is possible to derive an analogous version of (3.22) with boundary contributions along the seabed. With respect to this, in § B.2, we derive the relation

(3.23) \begin{eqnarray} \frac {\partial (\boldsymbol {R}_F - \boldsymbol {R}_B)}{\partial t} &=& \overline {\nabla } \left ({\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} - {\int _{-h}^{\eta }} \theta _3 \, {\textrm {d}z} \right ) + {\overline {\boldsymbol {u}}}^{\perp }_F \left ( \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F \right ) \nonumber \\ &&{}+ {\overline {\boldsymbol {u}}}^{\perp }_B \left ( \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B \right ) + \left [\overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F \, \overline {\nabla } \eta \right ] - \left [\overline {\boldsymbol {\theta }}_B - \left ( \theta _3 \right )_B \overline {\nabla } h\right ] . \end{eqnarray}

The above equation allows one to obtain a version of (3.22) with boundary conditions at the seabed. Indeed, subtracting (3.23) from (3.22), we immediately obtain

(3.24) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right ) &=& \overline {\nabla } \left ( w_F \, \frac {\partial \eta }{\partial t} - P_F - g \eta - {\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} + {\int _{-h}^{\eta }} \theta _3 \, {\textrm {d}z} \right ) \nonumber \\ && {}- {\overline {\boldsymbol {u}}}^{\perp }_B \left ( \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B \right ) + \left [\overline {\boldsymbol {\theta }}_B - \left ( \theta _3 \right )_B \, \overline {\nabla } h \right ] . \end{eqnarray}

Hereinafter, (3.22) and (3.24) will be the reference equations of the models studied in the present work. To recast them in closed forms, an expression for the flux $( w_F\, \partial \eta /\partial t - P_F)$ is finally needed. This is also important to avoid the presence of a time derivative inside the fluxes, since this is often a source of instability for numerical schemes. This alternative form is derived in Appendix C through the use of the kinematic condition at the free surface. In particular, we have the identity

(3.25) \begin{eqnarray} w_F \, \frac {\partial \eta }{\partial t} - P_F = \frac {N_F^2 \, w_F^2}{2} - \frac {1}{2} \left \|\frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right \|^2 - p_F . \end{eqnarray}

We recall that $w_F$ is computed using (2.4) after (3.5) for $\Upsilon$ is solved, $\boldsymbol {R}_F$ is provided by an external model for vorticity, and $p_F$ is obtained by introducing a proper closure for viscous/turbulent terms.

Before proceeding to the study of the derived models, a final remark is addressed. Equations (3.22) and (3.24) have to be compatible with (2.15) and (2.16) (which are only based on a decomposition of the velocity and vorticity fields and do not depend on any specific evolution equation). This implies that the following relations have to be satisfied:

(3.26) \begin{align} \frac {\partial (\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)}{\partial t} &= \overline {\nabla } \times \left [{\overline {\boldsymbol {u}}}^{\perp }_F \left ( \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F \right ) \right ] + \overline {\nabla } \times \left [\overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F\, \overline {\nabla } \eta \right ], \end{align}

(3.27) \begin{align} \frac {\partial (\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)}{\partial t} &= \overline {\nabla } \times \left [{\overline {\boldsymbol {u}}}^{\perp }_B \left ( \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B \right ) \right ] - \overline {\nabla } \times \left [\overline {\boldsymbol {\theta }}_B - \left ( \theta _3 \right )_B\, \overline {\nabla } h \right ] . \end{align}

The validity of these relations is proved in § B.3. Equations (3.26) and (3.27) highlight an important mechanism at the basis of the models described in the following sections. In particular, they show that the viscous effects at the interfaces (i.e. the terms depending on $\boldsymbol {\theta }$ ) act as sources for the generation of vorticity in the normal directions. From a different perspective, this means that the boundary conditions $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ are not completely ‘free’, but, except for their initial values, they evolve according to (3.26) and (3.27).

3.3.1. Circulation of the generalized velocity

As explained above, (3.26)–(3.27) imply the compatibility of (3.22)–(3.24) with (2.15)–(2.16). In particular, the latter allow for the derivation of interesting circulation theorems, which are a direct application of the Stokes’ theorem on the horizontal plane. Let us consider a domain $\mathcal {D} \subset \mathbb {R}^2$ with a regular boundary $\partial \mathcal {D}$ . Integrating (2.15) and (2.16) over $\mathcal {D}$ and applying the Stokes theorem, we immediately obtain

(3.28) \begin{align} \oint _{\partial \mathcal {D}}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right )\cdot \textrm{d}\ell &= \int _{\mathcal {D}} \left ( \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F \right )\, \textrm{d}V , \end{align}

(3.29) \begin{align} \oint _{\partial \mathcal {D}}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right ) \cdot \textrm{d}\ell &= - \int _{\mathcal {D}} \left ( \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B \right )\, \textrm{d}V . \end{align}

If the vorticity is zero and $\mathcal {D}$ is simply connected, then the circulations of $(\boldsymbol {M}/{d} + \boldsymbol {R}_F)$ and $(\boldsymbol {M}/{d} + \boldsymbol {R}_B)$ are identically null along any closed loop contained in such a domain, implying that they are gradient fields in $\mathcal {D}$ .

4. The model equations

The system of governing equations considered hereinafter is made of the continuity equation (3.4), one momentum equation among (3.22) and (3.24), and finally the Poisson equation for $\Upsilon$ . After proper models for the vorticity field and closures for the turbulent terms are provided, the system is in closed form, since the depth-averaged equations, namely (3.4), (3.22) and (3.24), provide the boundary conditions and forcing term for the Poisson equation, and, in turn, the solution of $\Upsilon$ allows us to evaluate the fluxes in the depth-averaged equations.

In the present section, we further inspect (3.22) and (3.24). In particular, we bring to light the existence of velocity potentials for both irrotational and rotational flows. We again highlight that the rotational and turbulent parts of the system are assumed to be known (e.g. through proper models for wave breaking and/or for friction along the seabed).

We also point out that (3.22) and (3.24) are obtained from the Navier–Stokes equations without any approximation. Hence the existence of velocity potentials for both irrotational and rotational flows is not a peculiarity of the present models, rather a general property of the wave propagation phenomenon. We observe, however, that in many Boussinesq and non-hydrostatic schemes, the evidence of this feature may be hidden or compromised by the approximations introduced during their derivation.

4.1. Irrotational fluid

Let us assume that the fluid is irrotational, i.e. $\boldsymbol {\omega } = 0$ during the evolution. Such a hypothesis also implies that no vorticity has to be diffused inside the fluid domain from the boundaries (both free surface and seabed). This can be attained only by assuming the fluid to be inviscid as a consequence of (3.26) and (3.27). These hypotheses lead to a substantial simplification of (3.22) and (3.24). In particular, they both coincide with the following simple expression

(4.1) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} \right ) = \overline {\nabla } \left (\frac {N_F^2 \, w_F^2}{2} - \frac {1}{2} \left \|\frac {\boldsymbol {M}}{d} \right \|^2 - g \eta \right ) \end{eqnarray}

that is obtained through the use of (3.25). The above expression implies that there exists a potential $\Phi (t,x,y)$ such that $\boldsymbol {M}/d = \overline {\nabla } \Phi$ . As shown in Antuono et al. (Reference Antuono, Valenza, Lugni and Colicchio2019), $\Phi = \phi _F$ , where $\phi (t,x,y,z)$ is the velocity potential, and (4.1) corresponds to a rearrangement of the Bernoulli equation at the free surface. From a different perspective, $\boldsymbol {M}/d$ is the velocity field obtained from the Zakharov potential $\phi _F$ (see Zakharov Reference Zakharov1968). In particular, (4.1) becomes

(4.2) \begin{eqnarray} \frac {\partial \phi _F}{\partial t} + \frac {\| \overline {\nabla } \phi _F \|^2}{2} - \frac {N_F^2 \, w_F^2}{2} + g \eta = C(t) , \end{eqnarray}

where $C(t)$ is a function of time only. Equation (4.2) is also used in the context of fully nonlinear potential flow models, as described, for example, in Yates & Benoit (Reference Yates and Benoit2015), Mohanlal et al. (Reference Mohanlal, Harris, Yates and Grilli2023) and Benoit, Zhang & Ma (Reference Benoit, Zhang and Ma2024). In analogy with these works, in Appendix C we prove the identity

(4.3) \begin{eqnarray} \frac {\partial \Upsilon }{\partial n} \bigg |_F = - N_F \, w_F . \end{eqnarray}

This means that the solution of (4.2) along with the Poisson equation for $\Upsilon$ can be regarded as a sort of Dirichlet-to-Neumann problem (see e.g. Lannes Reference Lannes2013).

Aside from the above considerations, we highlight that the Bernoulli equation at the free surface exhibits a substantial difference in comparison to (4.1) and (4.2). In the former case, the equation is evaluated along a moving interface (namely, the free surface), whereas in the latter case, the equations are formulated on a two-dimensional plane. This implies that the numerical implementation of (4.1) and (4.2) is considerably simpler.

4.2. Rotational fluid

If the fluid is rotational, then (3.22) and (3.24) represent different reformulations of the momentum equation of the Navier–Stokes system. In particular, (3.22) highlights the boundary contributions at the free surface coming from the vortical and viscous terms, whereas (3.24) shows those along the seabed. These boundary terms suggest which are the basic conditions that have to be specified along each interface, namely

(4.4) \begin{eqnarray} \boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F , \; \left [\overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F \, \overline {\nabla } \eta \right ] \quad \text {and} \quad \boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B , \; \left [\overline {\boldsymbol {\theta }}_B - \left ( \theta _3 \right )_B \, \overline {\nabla } h \right ] , \end{eqnarray}

according to the specific formulation at hand. As shown in the next subsubsection, there is no preferable choice, and each formulation has its own strong points and weaknesses. In any case, some care has to be paid when handling the pairs of conditions in (4.4). In fact, they are not independent but are related to each other through (3.26) and (3.27), respectively. An explanatory example is the problem of wave propagation over a frictional seabed. In that case, a proper closure has to be provided for the viscous boundary term $[\overline {\boldsymbol {\theta }}_B - ( \theta _3 )_B\, \overline {\nabla } h]$ , while the vortical boundary term $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ has to be obtained from (3.27). A simpler but interesting situation is described in the next section.

4.2.1. The rotational potentials

Let us assume that the boundary conditions $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ are null during the flow evolution, whereas the vorticity is generally non-null inside the fluid domain. Consistently with (3.26) and (3.27), and assuming that the two-dimensional domain is simply connected, there exist two potential functions $\chi$ and $\varsigma$ such that

(4.5) \begin{eqnarray} \overline {\boldsymbol {\theta }}_F + \left ( \theta _3 \right )_F\, \overline {\nabla } \eta = \overline {\nabla } \chi, \quad \overline {\boldsymbol {\theta }}_B - \left ( \theta _3 \right )_B\, \overline {\nabla } h = \overline {\nabla } \varsigma . \end{eqnarray}

Then (3.22) simplifies as

(4.6) \begin{eqnarray} \frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right ) = \overline {\nabla } \left (\frac {N_F^2 \, w_F^2}{2} - \frac {1}{2} \left \|\frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right \|^2 - g \eta - p_F + \chi \right ) , \end{eqnarray}

while (3.24) becomes

(4.7) \begin{eqnarray} &&\frac {\partial }{\partial t}\left ( \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right ) \nonumber \\ &&{}\quad = \overline {\nabla } \left (\frac {N_F^2 \, w_F^2}{2} - \frac {1}{2} \left \|\frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B \right \|^2 - g \eta - p_F - {\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} + {\int _{-h}^{\eta }} \theta _3 + \varsigma \right ) . \end{eqnarray}

The above results confirm that both $(\boldsymbol {M}/d + \boldsymbol {R}_F)$ and $(\boldsymbol {M}/d + \boldsymbol {R}_B)$ are in the form of a gradient field. Accordingly, there exist two potential functions $\Xi$ and $\Pi$ such that

(4.8) \begin{eqnarray} \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F = \overline {\nabla } \Xi , \quad \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_B = \overline {\nabla } \Pi . \end{eqnarray}

These satisfy the equations

(4.9) \begin{align} \frac {\partial \Xi }{\partial t} + \frac {\| \overline {\nabla } \Xi \|^2 }{2} - \frac {N_F^2 \, w_F^2}{2} + g \eta + p_F - \chi &= C_F(t) , \end{align}

(4.10) \begin{align} \frac {\partial \Pi }{\partial t} + \frac {\| \overline {\nabla } \Pi \|^2 }{2} - \frac {N_F^2 \, w_F^2}{2} + g \eta + p_F + {\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} - {\int _{-h}^{\eta }} \theta _3 - \varsigma &= C_B(t) , \end{align}

where $C_F$ and $C_B$ are functions of time only. Since the fluid is rotational, $\Xi$ and $\Pi$ are called rotational potentials hereinafter. Incidentally, we highlight that the potentials $\chi$ and $\varsigma$ are generally unknown at this level, unless related Poisson equations are solved or proper closures are provided.

Before proceeding, it is important to address some observations about the hypotheses introduced above. The condition $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F) = 0$ appears to be suited for long-crested waves and, more specifically, for spilling breaking events. Indeed, the particular geometry of the wave front is expected to approximately match the above condition. Incidentally, we observe that the wave-breaking model described in Antuono et al. (Reference Antuono, Lucarelli, Bardazzi and Brocchini2022) is compatible with it. In that work, good agreement with the experimental measurements was observed for spilling breaking waves that were essentially two-dimensional. Conversely, further studies should be addressed to inspect its accuracy in the presence of short-crested seas. Condition $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B) = 0$ seems a natural condition when the friction along the seabed is neglected. For both the cases described above, we should in theory specify the forms of $\chi$ and $\varsigma$ . Anyway, it seems a reasonable choice to assume them to be constant, so that the expressions in (4.5) are identically null.

5. The conservative-variable model

The systems derived up to here can be included in the family of primitive-variable models, i.e. those systems that are expressed in terms of the depth-averaged velocity field. In the specific case of systems (3.22) and (3.24), the primitive variables are the generalized velocity fields $(\boldsymbol {M}/d+\boldsymbol {R}_F)$ and $(\boldsymbol {M}/d+\boldsymbol {R}_B)$ . As shown in the present section, it is also possible to derive models that are in the form of conservation laws. These imply the conservation of certain global quantities (in the present case, the generalized mass flux) , therefore they are called conservative-variable models.

The derivation of the latter systems is rather interesting, since it requires the rearrangement of some vortical contributions that, as a consequence of the procedure, become implicit in the scheme. Despite the systems in conservative and primitive variables still being equivalent at the continuum, the rearrangement of the vortical terms may lead to non-negligible differences after numerical implementation.

As anticipated above, to derive the conservative model, we first need to isolate some vortical and viscous term from the expression (3.14) of $G_1$ . Specifically, we decompose it as

(5.1) \begin{eqnarray} G_1 = \widetilde {G}_1 - {\int _{-h}^{\eta }} q_1 \, {\textrm {d}z} + {\int _{-h}^{\eta }} \theta _1 \, {\textrm {d}z}, \end{eqnarray}

where $\widetilde {G}_1$ contains the remaining terms in (3.14). As shown in Appendix D, we can express the vortical term as

(5.2) \begin{eqnarray} {\int _{-h}^{\eta }} q_1 \, {\textrm {d}z} &=& \frac {\partial }{\partial x}\left ({\int _{-h}^{\eta }} \frac {u^2-v^2-w^2}{2} \, {\textrm {d}z} \right ) + \frac {\partial }{\partial y} \left ({\int _{-h}^{\eta }} u v \, {\textrm {d}z} \right ) \nonumber \\ && {}+u_F \, \frac {\partial \eta }{\partial t} + u_B \, \frac {\partial h}{\partial t} + \frac {\| {\boldsymbol {u}}_F \|^2}{2} \, \frac {\partial \eta }{\partial x} + \frac {\| {\boldsymbol {u}}_B \|^2}{2} \, \frac {\partial h}{\partial x} , \end{eqnarray}

while a straightforward application of the Leibniz rule leads to

(5.3) \begin{eqnarray} {\int _{-h}^{\eta }} \theta _1 \, {\textrm {d}z} &=& \frac {\partial }{\partial x}\left ({\int _{-h}^{\eta }} \mathbb {V}_{11} \, {\textrm {d}z} \right ) + \frac {\partial }{\partial y}\left ({\int _{-h}^{\eta }} \mathbb {V}_{12} \, {\textrm {d}z} \right ) - p_F \, \frac {\partial \eta }{\partial x} + \tau ^{(1)}_B , \end{eqnarray}

where

(5.4) \begin{eqnarray} \tau ^{(1)}_B = - \left ( \mathbb {V}_{11}\right )_B \, \frac {\partial h}{\partial x} - \left ( \mathbb {V}_{12}\right )_B \, \frac {\partial h}{\partial y} - \left ( \mathbb {V}_{13}\right )_B . \end{eqnarray}

Substituting the above expression in (3.15) and rearranging, we finally rewrite the equation in the form of a conservation law:

(5.5) \begin{eqnarray} \frac {\partial M_1}{\partial t} &+& \frac {\partial }{\partial x}\left [g\, \frac {d^2}{2} + \frac {d}{2} \left (\left \| \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right \|^2 - N_F^2 w_F^2 + U^2 - V^2 \right ) + {\int _{-h}^{\eta }} \frac { {\delta u}^2 - {\delta v}^2 - w^2 }{2}\, {\textrm {d}z}\!\! \right ] \nonumber \\ &+& \frac {\partial }{\partial y}\left ( d U V + {\int _{-h}^{\eta }} {\delta u} {\delta v} \, \textrm{d}z \right ) = \frac {\partial h}{\partial x} \left [ g d + \frac {1}{2} \left (\left \| \frac {\boldsymbol {M}}{d} + \boldsymbol {R}_F \right \|^2 - N_F^2 w_F^2 - \| {\boldsymbol {u}}_B \|^2 \right ) \right ] \nonumber \\ &+&\frac {\partial h}{\partial t} \left ( \frac {\partial \Upsilon _B }{\partial x} - w_B \, \frac {\partial h}{\partial x} \right ) + S_1 + T_1 ,\!\! \end{eqnarray}

where $S_1$ and $T_1$ contain the remaining vortical and turbulent terms, namely

(5.6) \begin{eqnarray} S_1 &=& - \frac {\partial }{\partial x}\left ({\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} q_3 \, {\textrm {d}}\zeta \right ) + \left ({\int _{-h}^{\eta }} q_3 \, {\textrm {d}z} \right ) \frac {\partial h}{\partial x}, \nonumber \\ T_1 &=& \frac {\partial }{\partial x}\left ({\int _{-h}^{\eta }} {\textrm {d}z} {\int _{z}^{\eta }} \theta _3 \,{\textrm {d}}\zeta - p_F d + {\int _{-h}^{\eta }} \mathbb {V}_{11} \, {\textrm {d}z} \right ) + \frac {\partial }{\partial y}\left ({\int _{-h}^{\eta }} \mathbb {V}_{12} \, {\textrm {d}z} \right ) \nonumber \\ && {}+\left (- {\int _{-h}^{\eta }} \theta _3 \, {\textrm {d}z} + p_F \right ) \frac {\partial h}{\partial x} + \tau ^{(1)}_B . \end{eqnarray}

The usual symmetry between vector components and spatial derivatives is applied to obtain the equation along the $y$ -direction, whereas the components along the vertical direction, namely $q_3$ and $\theta _3$ , remain unchanged. Again, we highlight that the viscous components are explicit in (5.5), whereas the vortical contributions from $q_1$ lead to the generation of the terms containing $\boldsymbol {U}$ and $\delta \boldsymbol {u}$ and thus become implicit in the conservative scheme.

Incidentally, we observe that (5.5) represents a simpler version of the models described in Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017, Reference Antuono, Lucarelli, Bardazzi and Brocchini2022).

5.1. The long-wave approximation

The long-wave approximation relies on the hypothesis that the characteristic wavelength of the propagation phenomenon is much larger than the reference water depth. Essentially, this corresponds to assuming that the higher-order derivatives in the equations give small contributions in comparison to the leading-order terms, and that therefore they can be neglected.

In the present formulation, this is equivalent to putting $\Upsilon \equiv 0$ , and dropping all the terms related to the vertical dynamics as well as any variation along the vertical direction. For the sake of simplicity, we also assume that the fluid is inviscid. Under these hypotheses, the formulation (5.5) reduces to the NSWEs in conservative form, i.e.

(5.7) \begin{eqnarray} \frac {\partial \left ( \boldsymbol {U} \right )}{\partial t} + \overline {\nabla } \cdot \left ({d} \boldsymbol {U} \otimes \boldsymbol {U} + \frac {g d^2}{2} \, \unicode {x1D7D9} \right ) = g d \, \overline {\nabla } h , \end{eqnarray}

while the formulation (3.22) gives

(5.8) \begin{eqnarray} \frac {\partial \boldsymbol {U}}{\partial t} + \overline {\nabla } \left (\frac {\| \boldsymbol {U} \|^2}{2} + g \eta \right ) - \boldsymbol {U}^\perp \left ( \frac {\partial V}{\partial x} - \frac {\partial U}{\partial y} \right ) = 0 , \end{eqnarray}

where the last factor on the left-hand side is the vorticity of the depth-averaged velocity field $\boldsymbol {U}$ . Note that such a term has been obtained by representing $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ through $\boldsymbol {U}$ , instead of assuming it to be an external forcing term. Expanding the gradient term and simplifying, we obtain the NSWEs in primitive variables, namely

(5.9) \begin{eqnarray} \frac {\partial \boldsymbol {U}}{\partial t} + \boldsymbol {U} \cdot \overline {\nabla } \boldsymbol {U} + g \, \overline {\nabla } \eta = 0. \end{eqnarray}

Equation (5.7) is recovered by multiplying (5.9) by $d$ and rearranging the terms through the use of the continuity equation (3.4).

The above findings prove that the NSWEs in two (horizontal) dimensions are intrinsically rotational. In particular, such equations include the vorticity in the vertical direction, whereas the other components are null. This also implies that they reduce to the irrotational formulation for one-dimensional propagation.

Incidentally, we highlight that the long-wave approximation is also at the basis of models like the Korteweg–de Vries and Serre–Green–Naghdi equations (see e.g. Lannes Reference Lannes2013). In the present work, however, these models are not considered, for the sake of brevity. In general, non-hydrostatic and Boussinesq schemes can be obtained from the present framework by introducing a proper modelling of the Poisson equation (3.5) – e.g. by employing a modal decomposition as in Raoult et al. (Reference Raoult, Benoit and Yates2016) or a multi-layer representation as in Ma et al. (Reference Ma, Shi and Kirby2012) and Pan et al. (Reference Pan, Kramer and Piggott2019) – or approximate solutions for $\Upsilon$ (e.g. polynomial expansions along the vertical direction as in Gobbi et al. Reference Gobbi, Kirby and Wei2000; Kim et al. Reference Kim, Lynett and Socolofsky2009).

6. Proof of concept

In this section, we consider some numerical applications to corroborate the main findings of the previous sections. The fluid is assumed to be inviscid and the turbulence effects negligible. The aim is to show that in absence of wave breaking and friction along the seabed, the wave motion is characterized by a null vorticity for the field $\boldsymbol {M}/d$ , while the vorticity of the depth-averaged field $\boldsymbol {U}$ is generally different from zero because of the interaction with the bathymetry.

Specifically, we consider the interaction of a non-breaking solitary wave with a submerged shoal placed over a constant-depth seabed. The shoal, which is inspired by the experiments described in the thesis of Chawla (Reference Chawla1995), consists of the top cut off portion of a sphere of radius 9.1m, whose centre is placed at $x=5\,\textrm{m}$ and $y=6.93\,\textrm{m}$ . The perimeter of the shoal is given by

(6.1) \begin{equation} (x-5)^2+(y-6.93)^2 = 2.57^2, \end{equation}

while the seabed bathymetry is given by

(6.2) \begin{equation} h = h_0 +8.73-\sqrt {82.81-(x-5)^2-(y-6.93)^2}, \end{equation}

where $h_0 = 0.45\,\textrm{m}$ is the constant water depth of the wave basin. The solitary wave height is $H = 0.0675\,\textrm{m}$ , corresponding to a nonlinearity parameter $\epsilon = H/h_0 = 0.15$ , and is generated through the ninth-order solution by Fenton (Reference Fenton1972). The numerical domain is $(x,y) \in [0,18]\times [0, 14.1]\,\textrm{m}^2$ , the inflow boundary is at $x=0$ , and open-boundary conditions are applied at $x=18\,\textrm{m}$ . Finally, wall conditions are imposed along the planes $y=0$ and $y=14.1\,\textrm{m}$ . We highlight that the present numerical domain is smaller than the one considered in the experimental campaign of Chawla (Reference Chawla1995), where regular wave trains were considered. In fact, differently from Chawla, we are not interested in refraction and wave–wave interaction phenomena, rather in the deformation and bending of the wave front. In the latter case, the propagation of a single solitary wave allows for a simpler and more straightforward analysis. We recall one more time that the wave is non-breaking and that no friction is considered along the seabed. Hence the field $\boldsymbol {M}/d$ is expected to be a gradient field, according to the theoretical findings of the previous sections.

A uniform Cartesian grid with $\Delta x = \Delta y = 0.03\,\textrm{m}$ is used in the horizontal direction, whereas a stretched grid with $\Delta z = 0.03\,\textrm{m}$ for $z \geq -0.12\,\textrm{m}$ increasing to $4\, \Delta z$ at the seabed is employed in the vertical direction. The still-water level is at $z = 0$ . The total number of solution points is approximately 3 700 000, and the simulation is run in the time interval $[0, 13]\,\textrm{s}$ .

Different numerical formulations are implemented in order to better highlight the generality of some results. For the schemes described in § 4.1, we consider the primitive-variable model based on the use of (4.1) and the potential-like model given by (4.2). In the latter case, (4.2) is integrated in time to obtain the potential $\phi _F$ , and a fourth-order central finite difference is used to compute its gradient (namely, $\boldsymbol {M}/d$ ). Conversely, (4.1) is integrated through a MUSCL–Hancock scheme with an approximate HLL Riemann solver. Due to the similarity, the numerical implementation complies with that used for the scheme described in Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017, Reference Antuono, Valenza, Lugni and Colicchio2019). The same approach is used to solve the conservative-variable model described in § 5. Furthermore, the scheme of Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017, Reference Antuono, Valenza, Lugni and Colicchio2019) is used as an additional check of the results. In all the above cases, the Poisson equation for $\Upsilon$ is discretized through second-order finite difference, and solved by using the library pARMS of parallel solvers for distributed sparse linear systems of equations described in Li, Saad & Sosonkina (Reference Li, Saad and Sosonkina2003). For the time integration, a fourth-order Adams–Bashforth–Moulton predictor–corrector scheme is adopted for all the models. It is worth noting that all the schemes that are in the form of conservation laws – namely, the primitive-variable model in (4.1), the conservative-variable model of § 5, and the scheme of Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017, Reference Antuono, Valenza, Lugni and Colicchio2019) – give similar results, as a consequence of the analogous structure and numerical implementation. For this reason, only the results of the model of § 5 are shown.

Figure 1 displays a three-dimensional view of the free-surface evolution at four selected time instants, namely $t = 6.14, 6.61, 7.07, 7.51\,\textrm{s}$ , obtained by using the potential-like model given by (4.2). As the solitary wave propagates, the central part of the wave front starts bending because of the slower velocity close to the shoal (top left panel), then its height increases as a consequence of shoaling and nonlinear effects (top right and bottom left panels). After the solitary wave overtakes the shoal, the wave height slightly decreases, and the front bending reduces (bottom right panel).

The above dynamics is proposed again in a two-dimensional view in figure 2. In this case, the results of the potential-like solver (contour fields) are compared with the outputs obtained through the conservative-variable model described in § 5 (contour levels in thick black lines). It is evident that both models predict practically the same evolution for $t = 6.14, 6.61\,\textrm{s}$ (top panels), while some negligible discrepancies arise later (bottom panels). The above dynamics has been also checked through the model described in Antuono et al. (Reference Antuono, Colicchio, Lugni, Greco and Brocchini2017, Reference Antuono, Valenza, Lugni and Colicchio2019), which essentially confirmed the same evolution.

Figure 1.

Snapshots of the free-surface evolution during the interaction of the solitary wave with the submerged shoal as computed by using (4.2) for $\phi _F$ .

Figure 2.

Snapshots of the free-surface evolution computed by using (4.2) for $\phi _F$ (contour fields) and the conservative-variable model described in § 5 (contour levels with thick black lines). The thin solid lines represent the contour levels of the seabed $h$ .

For the potential-like solver, figure 3 shows the vorticity fields computed from the depth-averaged velocity $\boldsymbol {U}$ (left column) and the generalized velocity $\boldsymbol {M}/d$ (right column). These are obtained through fourth-order central finite differences. Consistently with the theoretical findings of § 4.1, the field $\overline {\nabla } \times (\boldsymbol {M}/d)$ is zero up to the machine precision, while $\overline {\nabla } \times \boldsymbol {U}$ is sensibly different from zero in the regions where the wave front is characterized by the largest height and curvature. In particular, the magnitude of $\overline {\nabla } \times \boldsymbol {U}$ increases as the height and curvature increase (see panels at $t=6.14, 6.61\,\textrm{s}$ ), then decreases when the front bending reduces (see panels at $t=7.07, 7.51\,\textrm{s}$ ).

Figure 3.

Evolution of $\overline {\nabla } \times \boldsymbol {U}$ (left column) and $\overline {\nabla } \times (\boldsymbol {M}/d)$ (right column) at different time instants as computed by using (4.2) for $\phi _F$ . The thin solid lines represent the contour levels of the seabed $h$ .

The above analysis is repeated in figure 4 by using the outputs of the conservative-variable model of § 5. In this case, the overall dynamics is essentially the same, but both the fields $\overline {\nabla } \times \boldsymbol {U}$ and $\overline {\nabla } \times (\boldsymbol {M}/d)$ appear to be more noisy, especially close to the shoal boundary. This phenomenon is caused by the explicit presence of $\overline {\nabla } h$ in the source term. Furthermore, even if the magnitude of $\overline {\nabla } \times (\boldsymbol {M}/d)$ is smaller than $\overline {\nabla } \times \boldsymbol {U}$ , some spurious noise is still present inside the domain. This is as a consequence of the fact that $\boldsymbol {M}/d$ does not appear explicitly in the form of a gradient field, therefore its numerical computation is unavoidably affected by some inaccuracy. In addition to this issue, the use of the MUSCL–Hancock scheme implies that the spatial derivatives of the momentum fluxes are computed by using a variable reconstruction at the cell sides, and in turn, this implies that the post-processing application of a central finite difference to compute the curl of a gradient field generally leads to a non-null output.

Figure 4.

Evolution of $\overline {\nabla } \times \boldsymbol {U}$ (left column) and $\overline {\nabla } \times (\boldsymbol {M}/d)$ (right column) at different time instants as computed by using the conservative-variable model described in § 5. The thin solid lines represent the contour levels of the seabed $h$ .

Concluding, the above results confirm the theoretical findings of the previous sections, but at the same time suggest that a careful numerical implementation has to be considered in order to accurately detect the gradient components of the two-dimensional velocity fields.