A theory of depth averaging in models for coastal dynamics

Abstract

The present work proposes a general analysis of those models for gravity wave propagation that partially or totally rely on an average procedure over the water depth. The aim is the identification of the intrinsic physical quantities that characterize the wave dynamics, going beyond the usual definition of depth-averaged velocity. In particular, the proposed approach is based on the decomposition of the depth-averaged fields in their gradient- and divergence-free components. This naturally leads to the definition of a generalized velocity field that includes part of the dispersive contributions of the wave dynamics, and to the detection of the intrinsic boundary conditions along the free surface and the seabed. The analysis also proves the existence of generalized velocity potentials that under particular circumstances can include rotational contributions.

1. Introduction

The dynamics of gravity waves is a classic topic in fluid dynamics that captivated several generations of physicists, mathematicians and engineers. Its allure does not lie exclusively in the theoretical and physical aspects of wave propagation, but is also in the important practical applications, some directly related to human activities in the coastal and ocean regions. Different approaches can be employed for the theoretical modelling of this phenomenon according to the specific assumptions that are adopted. A complete description requires the implementation of the Navier–Stokes equations along with proper closures for turbulent terms (as in large-eddy simulations or Reynolds-averaged Navier–Stokes). In any case, if one assumes that turbulence plays a minor role, then it is possible to disregard the presence of boundary layers and dissipation, and rely on the use of the Euler equation. Moreover, if the vorticity in the fluid bulk is absent, and no singularities are present along the free surface and the seabed (namely, wave breaking and sharp seabed profiles), then the potential flow theory holds true and one can describe the wave dynamics through the velocity potential and the Bernoulli equation. All the approaches mentioned above allow for an accurate description of the wave dynamics, even though their applicability is limited to narrow sea regions when they are employed in numerical simulations. This issue is due to the fact that wave propagation generally affects large areas, whereas the above models depend on local quantities and are therefore more suited for the description of small-scale phenomena.

These aspects motivated the development of depth-averaged models, i.e. those models based on the assumption that the horizontal dynamics is predominant in the dynamics along the vertical direction and relying on the use of an average procedure over the water depth. These models can be obtained by applying the average procedure to each of the small-scale equations described above (i.e. Navier–Stokes equations, Euler equation, potential flow theory) along with proper theoretical assumptions. The simplest depth-averaged model is represented by the nonlinear shallow-water equations (NSWEs), which rely on the approximation of the pressure field through its hydrostatic component. Despite their simplicity and the availability of a wide collection of theoretical and numerical works (Carrier & Greenspan 1958; Stoker 1992; Toro 2001; Brocchini et al. 2001; LeVeque 2002; Antuono 2010; Rybkin, Pelinovsky & Didenkulova 2014), the accuracy of these equations is limited to the long-wave dynamics, i.e. wave propagation characterized by $h_0/\lambda \ll 1$ , where $h_0$ and $\lambda$ are respectively the reference water depth and wavelength. A more accurate approximation of the pressure field (with the inclusion of nonlinear and dispersive effects) leads to the definition of the Boussinesq-type models, which for many years have been a very active field of research (Madsen & Schäffer 1998; Veeramony & Svendsen 2000; Gobbi, Kirby & Wei 2000; Briganti et al. 2004; Eskilsson & Sherwin 2006; Kim, Lynett & Socolofsky 2009; Bingham, Madsen & Furman 2009; Tonelli & Petti 2012; Judge et al. 2018; Duran & Richard 2020). In comparison with the NSWEs, these models can describe the wave dynamics up to the intermediate-water conditions (namely $h_0/\lambda \lt 1/2$ ). The main drawback of Boussinesq-type models lies in their derivation, since this requires proper closures and/or specific approximations to obtain closed formulations in terms of depth-averaged variables. This limits their accuracy and consequently their range of applicability. A thorough description of Boussinesq models can be found, for example, in Brocchini (2013) and Kirby (2016).

For the purpose of overcoming the above limitations, many works in the last decades have addressed the definition of the so-called non-hydrostatic schemes (Zijlema & Stelling 2008; Ma, Shi & Kirby 2012; Raoult, Benoit & Yates 2016; Pan, Kramer & Piggott 2019). The basic idea is to represent, totally or partially, the vertical dynamics in order to recover a higher accuracy in the pressure field estimation. This allows non-hydrostatic schemes to describe the wave propagation up to the deep-water regime at the price of an increase of the computational costs in comparison to Boussinesq-type equations. Some non-hydrostatic schemes are obtained from equations that belong to large-eddy simulations (LES) or Reynolds-averaged Navier–Stokes (RANS) approaches (Christensen & Deigaard 2001; Smit, Zijlema & Stelling 2013; Derakhti et al. 2016), and for this reason are expected to share in part some of the computational limitations of the small-scale equations described before. Alternatively, other schemes are obtained by coupling depth-averaging procedures with a model (approximated or exact) for the vertical dynamics (Zijlema & Stelling 2005; Yamazaki, Kowalik & Cheung 2009; Antuono & Brocchini 2013; Lu & Xie 2016; Cantero-Chinchilla, Castro-Orgaz & Khan 2018; Pan et al. 2019; Escalante et al. 2024). In this case, the representation of the depth-averaged terms is again a crucial matter for investigation.

In comparison with small-scale equations, the models that totally or partially adopt a depth-averaging procedure cannot rely on a deepened physical characterization such as, for example, the potential flow theory. In this context, the main concern is about the closure and approximation of those terms that cannot be expressed directly as functions of depth-averaged quantities. Such a process, though accurate, has to be applied with great care, since it may lead to the detriment of general physical properties in the models. Furthermore, many approaches are driven by the necessity of attaining numerical schemes characterized by limited computational costs, rather than a thorough physical analysis. In fact, there is not yet a general and comprehensive overview about the physical structure that such schemes have to possess.

The aim of the present work is therefore to provide a contribution on this topic, and derive a theory for those models that (partially or totally) employ depth-averaging procedures. This is done by identifying the intrinsic physical quantities that characterize the wave dynamics through the decomposition of the averaged quantities in their gradient- and divergence-free components. This allows for a clear identification of the rotational and irrotational terms, and likewise the local models. Furthermore, the adopted approach allows for a direct inspection of the necessary conditions along the interfaces, namely the free surface and the bottom seabed.

The paper is organized as follows. Section 2 introduces some basic concepts that are used in the derivation that follows, while in § 3, we derive the general depth-averaged equations starting from the Navier–Stokes equations. Hence an extensive analysis of the governing equations is provided in Sections 4 and 5, while a proof of concept of the main theoretical findings is provided in § 6. Finally, an overall discussion of the results is proposed in § 7.

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. (2017), 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. (1995), Erduran, Ilic & Kutija (2005) and Tonelli & Petti (2012), 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 1993; Chen et al. 2000; Kennedy et al. 2000), 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. (2022), 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. (2017). 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 1993; Kim et al. 2009; Donahue et al. 2015; Kazakova & Richard 2019) and/or to the Laplace equation for the velocity potential $\phi$ (e.g. Wei et al. 1995; Gobbi et al. 2000; Raoult et al. 2016). 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. (2019) 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. (2019), $\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 1968). 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 (2015), Mohanlal et al. (2023) and Benoit, Zhang & Ma (2024). 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 2013).

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. (2022) 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. (2017, 2022).

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 2013). 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. (2016) or a multi-layer representation as in Ma et al. (2012) and Pan et al. (2019) – or approximate solutions for $\Upsilon$ (e.g. polynomial expansions along the vertical direction as in Gobbi et al. 2000; Kim et al. 2009).

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 (1995), 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 (1972). 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 (1995), 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. (2017, 2019). The same approach is used to solve the conservative-variable model described in § 5. Furthermore, the scheme of Antuono et al. (2017, 2019) 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 (2003). 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. (2017, 2019) – 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. (2017, 2019), 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.

7. Discussion and conclusions

The analysis proposed so far highlights some features that a model which partially or totally relies on a depth-average procedure has to exhibit, but at the same time, it clears the way to further questions and lines of research. In the following discussion, we address these along with a brief summary of the results obtained up to here.

In § 3.3, we showed that starting from the Navier–Stokes equations, it is possible to derive two equivalent depth-averaged formulations with fluxes that are in the form of gradient fields, with forcing terms containing contributions along the free surface or along the seabed. These formulations have been derived under general assumptions and without any approximation, therefore can be regarded as a rearrangement of the Navier–Stokes equations. The analysis performed on these depth-averaged models highlights a number of interesting aspects that are important for the physical characterization of the wave propagation phenomenon. In particular, we showed that the classic depth-averaged velocity is not the most significant field for the description of the wave dynamics. The depth-averaged system, in fact, is naturally well represented through a generalized velocity field that partially includes the dispersive effects of wave propagation. This also allowed us to highlight some intriguing physical features. Specifically, we identified the conditions under which the existence of generalized velocity potentials can be established. Differently from the classical potential theory, the generalized potentials can also include rotational contributions, and stem from a more general decomposition of the depth-averaged fields in their gradient- and divergence-free components.

As a proof of concept, the models described in §§ 4.1 and 5 have been implemented numerically and applied to simulate the evolution of a non-breaking solitary wave interacting with a frictionless uneven bottom. The numerical outputs proved to be consistent with the theoretical findings, even though they suggest that the identification of the gradient fields is a delicate and insidious matter. The diffculty is essentially related to two interconnected aspects, namely the specific structure of the model under consideration and its numerical implementation. In particular, the use of (4.2) for $\phi _F$ , or alternatively of the primitive-variable formulation (4.1), seems the most reliable approach, even though in the latter case, numerical approximations may jeopardize the results. The implementation of ad hoc numerical techniques is, however, expected to overcome this issue.

The above considerations are useful for the definition of future schemes for coastal dynamics. In fact, the structures of the models described in the present work can be used as general frameworks for the derivation of non-hydrostatic or Boussinesq schemes that prescribe the presence of generalized fields (namely, $\boldsymbol {M}$ and $\boldsymbol {M}/d$ ) and gradient fields (under the conditions reported in §§ 4.1 and 4.2.1). As a consequence of the results shown in §§ 4.1 and 5.1, the generalized fields are expected to be good candidates for the description the dynamics in both deep- and shallow-water conditions. Indeed, being closely related to the potential $\phi _F$ , they remain significant in deep-water conditions, whereas $\boldsymbol {U}$ loses its physical relevance, since the depth-average procedure is applied for a large extent to regions where the local fluid motion is weak. Vice versa, the generalized fields naturally tend to depth-averaged quantities when the long-wave approximation is valid, as shown in § 5.1. For the modelling of gradient fields, this seems an important aspect if one aims at a correct description of the vorticity. Indeed, an inaccurate representation of these terms may lead to a non-physical generation of vorticity at the interfaces as a consequence of (2.15) and (2.16), and to its subsequent diffusion inside the fluid bulk.

An essential subject that still needs to be inspected is the match of the formulations described in (3.22) and (3.24) with models for the dynamics of the vorticity field. The latter have to include a description of the breaking process and of the consequent generation of turbulence. In this case, a study on the essential boundary conditions for the vorticity field and their link to the boundary conditions of the depth-averaged models would be beneficial for the understanding of the wave dynamics phenomenon.

A further important point to address in future is about the numerical implementation of the two alternative formulations represented in (3.22) and (3.24) for breaking waves and/or frictional seabeds, since they may exhibit different pros and cons according to the problem at hand.

Finally, we observe that many aspects described in the previous sections drastically simplify in two dimensions (i.e. for unidirectional wave propagation). In particular, the terms $(\boldsymbol {\omega }_F \cdot \boldsymbol {\tilde {n}}_F)$ and $(\boldsymbol {\omega }_B \cdot \boldsymbol {\tilde {n}}_B)$ are identically null in two dimensions, as well as (3.26) and (3.27). Conversely, the viscous terms $[\overline {\boldsymbol {\theta }}_F + ( \theta _3 )_F\, \overline {\nabla } \eta ]$ and $[\overline {\boldsymbol {\theta }}_B - ( \theta _3 )_B\, \overline {\nabla } h]$ are generally different from zero and need to be modelled. These aspects, which are not simply related to a higher geometrical complexity but include specific physical features, confirm that the three-dimensional case is intrinsically more rich than the two-dimensional one.