2. The 2-D long-wave PCSSE

The equations governing the irrotational flow of an incompressible, inviscid fluid with a free surface, in terms of the polar coordinates, are given by

(2.1) \begin{equation} \frac {1}{r} \frac {\partial }{\partial r} \left ( r \frac {\partial \psi }{\partial r} \right ) + \frac {1}{r^2} \frac {\partial ^2 \psi }{\partial \theta ^2} = 0, \quad-\alpha \lt \theta \lt 0, \end{equation}

(2.2) \begin{equation} \frac {1}{r} \frac {\partial \psi }{\partial \theta } +\frac {1}{\sigma } \frac {\partial ^2 \psi }{\partial r^2}= 0,\quad \theta = 0, \end{equation}

(2.3) \begin{equation} \psi = 0, \quad \theta = -\alpha , \end{equation}

where $r$ and $\theta$ are the radial distance and angle, respectively, $\theta =0$ represents the mean free surface, $\theta =-\alpha$ the bottom location and $r=0$ the shore, as depicted in figure 1. This formulation refers to a time-harmonic flow such that the time-harmonic stream function is defined by $\textrm{Re}[\psi \textrm{e}^{-\textrm{i} \omega t}]$ , where $\omega$ is the angular frequency of the incoming wave and $\sigma = {\omega ^2}/{g}$ , $g$ being the gravitational acceleration. Function $\psi = \psi (r,\theta )$ is a complex-valued function.

Schwartz et al. (Reference Schwartz, Oron and Agnon2023) derived the PCSSE using the following normalised, $\theta$ -dependent profile:

(2.4) \begin{equation} F_p(\theta ) = \frac {\sinh(\gamma (\theta +\alpha ))}{\sinh(\gamma \alpha )}, \end{equation}

Figure 1. Polar coordinates: problem description.

where $\gamma$ is a constant, equivalent to the wavenumber in the Cartesian formulation. The vertical analytical integration of the Cartesian-equivalent profile produced complex terms, necessitating the numerical calculation of the coefficients. In the following, we consider the long-wave limit of this profile, which significantly simplifies the coefficients and enables the use of analytical expressions based on this approximation. This is achieved by taking the leading-order term of the polar profile for a long-wave approximation ( $\gamma \theta , \gamma \alpha \ll 1$ ) (2.4):

(2.5) \begin{equation} F_p(\theta ) \approx 1 + \frac {\theta }{\alpha } + O \left ((\gamma \alpha )^3\right ). \end{equation}

The equivalent Cartesian terms are derived by substituting the following relations (figure 2):

(2.6) \begin{equation} \theta = \arctan\left (\frac {z h'}{h}\right ), \quad \alpha = \arctan(h'). \end{equation}

Here $h'$ is the slope of the bottom bathymetry and is defined in figure 1. The vertical Cartesian-equivalent profile, neglecting the higher-order terms, is given as

(2.7) \begin{equation} F_{pc}\left (h(x),h'(x),z\right ) = 1+\frac {\arctan\left (z h'/h\right )}{\arctan(h')} . \end{equation}

Figure 2. Polar–Cartesian relations.

The vertical integration of this simplified profile still presents some complications. Integrating $\arctan(n)$ with respect to a general parameter $n$ leads to expressions involving non-elementary functions. To avoid these terms and simplify the mathematical representation of the coefficients, we utilise an approximation of this function, as provided in Hastings (Reference Hastings1955). This approximation is represented as a polynomial of $n$ and has a maximum error of 0.6 per mil in the range $\mid n \mid \leq 1$ :

(2.8) \begin{equation} App_1{:}\quad \arctan(n) \approx c_1 n + c_3 n^3 + c_5 n^5 ,\end{equation}

with $c_1 = 0.995354, \hspace {2mm} c_3 = -0.288679$ and $c_5 = 0.079331$ . In the rest of the region, a similar approximation for the inverse parameter ${1}/{n}$ can be utilised, providing a rational function approximation:

(2.9) \begin{align} &\quad\!\! App_2{:}\quad \arctan(n) \approx \frac {Pi}{2} - \left (c_1 \frac {1}{n} + c_3 \frac {1}{n^3} + c_5 \frac {1}{n^5}\right ) \quad \text{in} \ n\gt 1 , \end{align}

(2.10) \begin{align} &App_3{:}\quad \arctan(n) \approx -\frac {Pi}{2} - \left (c_1 \frac {1}{n} + c_3 \frac {1}{n^3} + c_5 \frac {1}{n^5}\right ) \quad \text{in}\ n\lt -1 \end{align}

with the same maximum error.

Noting that $({z}/{h}) \leq 1$ , the term in the numerator of (2.7), $(zh'/h) $ , can be classified by dividing it into three regions:

1. (i) $\mid h' \mid \leq 1$ . The absolute value of the integrand is smaller than that along the integration range, and therefore: $\int _{-h}^{0} F \,\textrm{d}z = \int _{-h}^{0} App_1 \,\textrm{d}z$ .

2. (ii) $h'\gt 1$ . The integrand crosses two approximation regions, and therefore the integration is divided into two parts: $\int _{-h}^{0} F \,\textrm{d}z = \int _{-h}^{-h/h'} App_3 \,\textrm{d}z + \int _{-h/h'}^{0} App_1 \,\textrm{d}z$ .

3. (iii) $h'\lt -1$ . The approximation is also divided into two parts: $\int _{-h}^{0} F \,\textrm{d}z = \int _{-h}^{h/h'} App_2 \,\textrm{d}z + \int _{h/h'}^{0} App_1 \,\textrm{d}z$ .

Applying this approximation requires distinguishing between the different regions of the bottom slope along the flow. Once the region is identified, the appropriate integrand can be substituted. An additional option is to numerically integrate the vertical profile. In both options, there are no limitations on the steepness of the bottom slope, leading to a 2-D steep-slope equation, referred to as the LPCSSE.

3. The 3-D PCSSE

The challenge in deriving an equation for a 3-D geometry stems from the varying direction of the depth gradient at each location. In addition, the wavenumber components in the $x$ and $y$ directions are influenced by the depth gradient direction, necessitating their determination at each point. By taking the long-wave limit, the velocity profiles become independent of the wavenumber, helping to overcome this challenge. The remaining issue is addressed by defining a local axis system where the xaxis is aligned with the depth gradient direction at each point and the $y$ axis is perpendicular to the $x$ axis. This results in a locally horizontal depth in the direction of $y$ .

The time-harmonic vector stream function is defined as

(3.1) \begin{equation} \boldsymbol{\Psi }(x,y,z) = (f_1 \psi _1(x,y), f_2 \psi _2(x,y)). \end{equation}

The vertical profiles of the vector stream function are defined by

(3.2) \begin{align} &f_1(h,h',z) = \frac {\sinh\left (\gamma (\arctan(z h'/h) + \arctan(h')\right )}{\sinh \left (\gamma \arctan(h')\right )}, \end{align}

(3.3) \begin{align} &\qquad\qquad\qquad f_2(h,z) = \frac {\sinh\left (k(z+h)\right )}{\sinh(k h)} , \end{align}

where $\gamma$ is the $x$ component and $k$ is the $y$ component of the wavenumber vector, $f_1$ is the Cartesian equivalent of the polar profile and $f_2$ corresponds to a horizontal bottom.

The velocities are calculated accordingly:

(3.4) \begin{align} &\mathbf{u}(x,y,z) = \left (u(x,y,z), v(x,y,z)\right ) = \left (f_{1_z} \psi _1, f_{2_z} \psi _2\right ), \end{align}

(3.5) \begin{align} &\quad w(x,y,z) = -\nabla \cdot \boldsymbol{\Psi } = -f_{1_{x}} \psi _1 - f_1 \psi _{1_{x}} - f_2 \psi _{2_{y}}, \end{align}

where $\nabla = (\partial _{x},\partial _{y})$ .

The free-surface location is derived from the linear kinematic condition:

(3.6) \begin{equation} \eta (x,y) = \frac {1}{\textrm{i}\omega } \nabla \cdot \varPsi _0(x,y) = \frac {1}{\textrm{i}\omega } \left (\psi _{1_{x}}(x,y) + \psi _{2_{y}}(x,y)\right ). \end{equation}

These terms are substituted into the time-averaged Lagrangian density:

(3.7) \begin{align} L &= \int _{-h}^{0} \left (|u|^2+|v|^2+|w|^2\right )\textrm{d}z - g |\eta |^2 \notag \\ & = d_1 |\psi _1|^2 + d_2 |\psi _2|^2 + c |\psi _1|^2 + a_1 |\psi _{1_{x}}|^2 + a_2 |\psi _{2_{y}}|^2 \notag \\ &\quad + 2 b_1 \textrm{Re}(\psi _1 \psi _{1_{x}}*) + 2 b_2 \textrm{Re}(\psi _1 \psi _{2_{y}}*) + 2 e \textrm{Re}(\psi _{1_{x}} \psi _{2_{y}}*), \end{align}

where the coefficients are defined as

(3.8) \begin{align} &a_1 = \int f_1^2 \,\textrm{d}z - \frac {1}{\sigma }, \quad a_2 = \int f_2^2 \,\textrm{d}z - \frac {1}{\sigma }, \quad b_1 = \int f_1 f_{1_{x}}\, \textrm{d}z, \quad b_2 = \int f_2 f_{1_{x}} \,\textrm{d}z, \end{align}

(3.9) \begin{align} &\qquad c = \int f_{1_{x}}^2 \,\textrm{d}z, \quad e = \int f_1 f_2\, \textrm{d}z - \frac {1}{\sigma }, \quad d_1 = \int f_{1_z}^2\, \textrm{d}z, \quad d_2 = \int f_{2_z}^2 \,\textrm{d}z. \end{align}

The variation with respect to $\psi _1$ results in

(3.10) \begin{align} \frac {\delta L}{\delta \psi _1} = \frac {\partial L}{\partial \psi _1} - \frac {\partial }{\partial x} \frac {\partial L}{\partial \psi _{1_{x}}} = 0, \nonumber \\ a_1 \psi _{1_{xx}} + e \psi _{2_{xy}} + a_{1_{x}} \psi _{1_{x}} + (e_{x} - b_2) \psi _{2_{y}} + (b_{1_{x}}-c -d_1) \psi _1 = 0. \end{align}

The variation with respect to $\psi _2$ provides the coupled equation

(3.11) \begin{align} \frac {\delta L}{\delta \psi _2} = \frac {\partial L}{\partial \psi _2} - \frac {\partial }{\partial y} \frac {\partial L}{\partial \psi _{2_{y}}} = 0, \nonumber \\ a_2 \psi _{2_{yy}} + e \psi _{1_{xy}} + b_2 \psi _{1_{y}} - d_2 \psi _2 = 0. \end{align}

Figure 3. The dynamic $x$ , $y$ and static $\tilde {x}$ , $\tilde {y}$ Cartesian axis systems. The $x$ axis is aligned with the direction of the local bottom gradient. The $z$ and $\tilde {z}$ axes are parallel to the drawing plane. Angle $\delta$ represents the counterclockwise angle between the two axis systems.

Equations (3.10) and (3.11) constitute a coupled system on the dynamic $x,y$ axes. The equations are now rotated to the static axis system for which $\tilde {x}$ is perpendicular and $\tilde {y}$ is parallel to the shoreline. Defining $\delta$ as the counterclockwise rotation angle between the two systems (figure 3), we can define

(3.12) \begin{align} \tilde {x} = x \cos {\delta } - y \sin {\delta }, \end{align}

(3.13) \begin{align} \tilde {y} = x \sin {\delta } + y \cos {\delta }. \end{align}

The relationships between the derivatives in the $(\tilde {x},\tilde {y})$ and $(x,y)$ coordinate systems are given by

(3.14) \begin{align}&\qquad\qquad\quad \psi _{{y}} = -\sin (\delta ) \tilde {\psi }_{\tilde {x}} +\cos (\delta ) \tilde {\psi }_{\tilde {y}}, \end{align}

(3.15) \begin{align} &\qquad \psi _{{yy}} = \sin ^2(\delta ) \tilde {\psi }_{{\tilde {x}\tilde {x}}} - \sin (2 \delta ) \tilde {\psi }_{{\tilde {x}\tilde {y}}} + \cos ^2(\delta ) \tilde {\psi }_{{\tilde {y}\tilde {y}}}, \end{align}

(3.16) \begin{align} &\qquad\qquad\qquad\psi _{{x}} = \cos (\delta ) \tilde {\psi }_{\tilde {x}} + \sin (\delta ) \tilde {\psi }_{\tilde {y}}, \end{align}

(3.17) \begin{align} &\qquad\psi _{{xx}} = \cos ^2(\delta ) \tilde {\psi }_{{\tilde {x}\tilde {x}}} + \sin (2 \delta ) \tilde {\psi }_{{\tilde {x}\tilde {y}}} + \sin ^2(\delta ) \tilde {\psi }_{{\tilde {y}\tilde {y}}}, \end{align}

(3.18) \begin{align} &\psi _{{xy}} = \frac {\sin (2 \delta )}{2} (\tilde {\psi }_{{\tilde {y}\tilde {y}}}-\tilde {\psi }_{{\tilde {x}\tilde {x}}}) + \left (\cos ^2(\delta )-\sin ^2(\delta ) \right ) \tilde {\psi }_{{\tilde {x}\tilde {y}}}. \end{align}

Additionally, the relationships between the coefficients are given by

(3.19) \begin{align} &b_1 = \cos (\delta ) \widetilde {bx}_1 + \sin (\delta ) \widetilde {by}_1, \quad b_2 = \cos (\delta ) \widetilde {bx}_2 + \sin (\delta ) \widetilde {by}_2, \end{align}

(3.20) \begin{align} &\qquad\qquad c = \cos ^2(\delta ) \widetilde {cx} + \sin (2\delta ) \widetilde {cxy} + \sin ^2(\delta ) \widetilde {cy}, \end{align}

where

(3.21) \begin{align} &\widetilde {bx}_1 = \int \tilde {f_1} \tilde {f}_{1_{\tilde {x}}}\, \textrm{d}z, \quad \widetilde {by}_1 = \int \tilde {f_1} \tilde {f}_{1_{\tilde {y}}}\, \textrm{d}z, \quad \widetilde {bx}_2 = \int \tilde {f_2} \tilde {f}_{1_{\tilde {x}}}\, \textrm{d}z, \quad \widetilde {by}_2 = \int \tilde {f_2} \tilde {f}_{1_{\tilde {y}}}\, \textrm{d}z, \end{align}

(3.22) \begin{align} &\qquad\qquad\widetilde {cx} = \int \tilde {f}_{1_{\tilde {x}}}^2\, \textrm{d}z, \quad \widetilde {cxy} = \int \tilde {f}_{1_{\tilde {x}}} \tilde {f}_{1_{\tilde {y}}}\, \textrm{d}z, \quad \widetilde {cy} = \int \tilde {f}_{1_{\tilde {y}}}^2\, \textrm{d}z. \end{align}

Equations (3.10) and (3.11) are now expressed in the static coordinate system. The transformed form of equation (3.10) is given by

(3.23) \begin{align} & \tilde {a}_1 \left (\cos ^2(\delta ) \tilde {\psi }_{1_{\tilde {x}\tilde {x}}} + \sin (2 \delta ) \tilde {\psi }_{1_{\tilde {x}\tilde {y}}} + \sin ^2(\delta ) \tilde {\psi }_{1_{\tilde {y}\tilde {y}}}\right ) \notag \\ &\quad + \tilde {e} \left (\frac {1}{2} \sin (2 \delta )(\tilde {\psi }_{2_{\tilde {y}\tilde {y}}}-\tilde {\psi }_{2_{\tilde {x}\tilde {x}}}) + (\cos ^2(\delta )-\sin ^2(\delta )) \tilde {\psi }_{2_{\tilde {x}\tilde {y}}}\right ) \notag \\ &\quad + \left (\cos (\delta ) \tilde {a}_{1_{\tilde {x}}} + \sin (\delta ) \tilde {a}_{1_{\tilde {y}}}\right ) \left (\cos (\delta ) \tilde {\psi }_{1_{\tilde {x}}} + \sin (\delta ) \tilde {\psi }_{1_{\tilde {y}}} \right ) \notag \\ & \quad+ \left (\cos (\delta ) \tilde {e}_{\tilde {x}} + \sin (\delta ) \tilde {e}_{\tilde {y}} - \cos (\delta ) \widetilde {bx}_2 - \sin (\delta ) \widetilde {by}_2\right ) \left (\cos (\delta ) \tilde {\psi }_{2_{\tilde {y}}} - \sin (\delta ) \tilde {\psi }_{2_{\tilde {x}}} \right ) \notag \\ &\quad + \left (\cos ^2(\delta ) \widetilde {bx}_{1_{\tilde {x}}} + \frac {1}{2} \sin (2 \delta ) \widetilde {by}_{1_{\tilde {x}}} + \frac {1}{2} \sin (2 \delta ) \widetilde {bx}_{1_{\tilde {y}}} + \sin ^2(\delta ) \widetilde {by}_{1_{\tilde {y}}} \right. \notag \\ &\quad -\left. \cos ^2(\delta ) \tilde {c}_{\tilde {x}} - \sin (2 \delta ) \tilde {c}_{\tilde {x} \tilde {y}} - \sin ^2(\delta ) \tilde {c}_{\tilde {y}} - \tilde {d}_1 \right )\tilde {\psi }_1 = 0. \end{align}

Similarly, the transformed form of equation (3.11) is given by

(3.24) \begin{align} & \tilde {a}_2 \left (\sin ^2(\delta ) \tilde {\psi }_{{2_{\tilde {x} \tilde {x}}}} - \sin (2 \delta ) \tilde {\psi }_{{2_{\tilde {x} \tilde {y}}}} + \cos ^2(\delta ) \tilde {\psi }_{{2_{\tilde {y} \tilde {y}}}}\right ) \notag \\ &\quad + \tilde {e} \left (\frac {1}{2} \sin (2 \delta )(\tilde {\psi }_{1_{\tilde {x} \tilde {y}}}-\tilde {\psi }_{1_{\tilde {x} \tilde {x}}}) + (\cos ^2(\delta )-\sin ^2(\delta )) \tilde {\psi }_{1_{\tilde {x} \tilde {y}}}\right ) \notag \\ &\quad + \left (\cos (\delta ) \widetilde {bx}_2 + \sin (\delta ) \widetilde {by}_2\right ) \left (\tilde {\psi }_{1_{\tilde {y}}} \cos (\delta ) - \tilde {\psi }_{1_{\tilde {x}}} \sin (\delta ) \right ) - \tilde {d}_2 \tilde {\psi }_2 = 0. \end{align}

Dividing both equations by $\cos ^2(\delta )$ provides

(3.25) \begin{align} & \tilde {a}_1 \left (\tilde {\psi }_{{1_{\tilde {x} \tilde {x}}}} + 2 \tan (\delta ) \tilde {\psi }_{{1_{\tilde {x} \tilde {y}}}} + \tan ^2(\delta ) \tilde {\psi }_{{1_{\tilde {y} \tilde {y}}}}\right ) \notag \\ &\quad + \tilde {e} \left (\tan (\delta )(\tilde {\psi }_{2_{\tilde {y} \tilde {y}}}-\tilde {\psi }_{2_{\tilde {x} \tilde {x}}}) + (1-\tan ^2(\delta )) \tilde {\psi }_{2_{\tilde {x} \tilde {y}}}\right ) \notag \\ &\quad + \left (\tilde {a}_{1_{\tilde {x}}} + \tan (\delta ) \tilde {a}_{1_{\tilde {y}}}\right ) \left (\tilde {\psi }_{1_{\tilde {x}}} + \tan (\delta ) \tilde {\psi }_{1_{\tilde {y}}} \right ) \notag \\ &\quad + \left (\tilde {e}_{\tilde {x}} + \tan (\delta ) \tilde {e}_{\tilde {y}} - \widetilde {bx}_2 - \tan (\delta ) \widetilde {by}_2\right ) \left (\tilde {\psi }_{2_{\tilde {y}}} - \tan (\delta ) \tilde {\psi }_{2_{\tilde {x}}} \right ) \notag \\ &\quad + \left (\widetilde {bx}_{1_{\tilde {x}}} + \tan (\delta ) \widetilde {by}_{1_{\tilde {x}}} + \tan (\delta ) \widetilde {bx}_{1_{\tilde {y}}} + \tan ^2(\delta ) \widetilde {by}_{1_{\tilde {y}}}\vphantom{\frac {\tilde {d}_1}{\tan ^2(\delta )}}\right. \notag \\ &\quad \left.- \tilde {c}_{\tilde {x}} - 2 \tan (\delta ) \tilde {c}_{\tilde {x} \tilde {y}} - \tan ^2(\delta ) \tilde {c}_{\tilde {y}} - \frac {\tilde {d}_1}{\tan ^2(\delta )} \right ) \tilde {\psi }_1 = 0, \end{align}

(3.26) \begin{align} & \tilde {a}_2 \left (\tan ^2(\delta ) \tilde {\psi }_{{2_{\tilde {x} \tilde {x}}}} - 2 \tan (\delta ) \tilde {\psi }_{{2_{\tilde {x} \tilde {y}}}} + \tilde {\psi }_{{2_{\tilde {y} \tilde {y}}}}\right ) \notag \\ &\quad + \tilde {e} \left (\tan (\delta )(\tilde {\psi }_{1_{\tilde {x} \tilde {y}}}-\tilde {\psi }_{1_{\tilde {x} \tilde {x}}}) + (1-\tan ^2(\delta )) \tilde {\psi }_{1_{\tilde {x} \tilde {y}}}\right ) \notag \\ &\quad + \left (\widetilde {bx}_2 + \tan (\delta ) \widetilde {by}_2\right ) \left (\tilde {\psi }_{1_{\tilde {y}}} - \tan (\delta ) \tilde {\psi }_{1_{\tilde {x}}} \right ) - \frac {\tilde {d}_2}{\cos ^2(\delta )} \tilde {\psi }_2 = 0. \end{align}

Equations (3.10) and (3.11) are to be solved in the dynamic $x$ , $y$ , $z$ coordinate system with the appropriate lateral boundary conditions. Following this, the system is rotated at each point by the rotation angle $\delta$ , as outlined in equations (3.12) and (3.13). Alternatively, equations (3.25) and (3.26) can be solved along with the rotated lateral boundary conditions.

4. The time-domain equation

Generally, mild-slope-type equations are time-harmonic, meaning they are defined in the frequency domain. Porter (Reference Porter2019) developed an ESWE by accounting for the next order in the wavelength. Following a similar procedure and using the PC vertical profiles, a MSE is formulated in the time domain. The 2-D case is addressed in § 4.1, while the 3-D case is covered in § 4.2.

4.1. The 2-D equation

The horizontal velocity is defined by a separation of variables:

(4.1) \begin{equation} u(x,z,t) = f_{1_z}(z,h(x),h'(x)) u_0(x,t), \end{equation}

where $f_1$ is the vertical profile as calculated in equation (3.2).

Integrating vertically provides the stream function

(4.2) \begin{equation} \psi (x,z,t) = f_1 u_0, \end{equation}

and taking the horizontal derivative provides the vertical velocity component

(4.3) \begin{equation} w(x,z,t) = -f_{1_x} u_0 + f_1 \zeta _t, \end{equation}

where the linear kinematic free-surface condition, $u_{0_x} = -\zeta _t$ , was substituted. Here, $\zeta (x,t) = \eta (x) \textrm{e}^{-\textrm{i} \omega t}$ .

It is noted that the bottom boundary condition, $w + h' u = 0 \ \text{on}\ z = -h$ , is satisfied under this representation. This is in addition to the continuity equation, which is inherently satisfied by the stream function definition.

Next, the horizontal component of the linearised momentum equation is utilised:

(4.4) \begin{equation} \rho w_t = -p_z - \rho g, \end{equation}

and by substituting the vertical velocity profile, we derive

(4.5) \begin{equation} p_z = \rho \left (f_{1_x} u_{0_t} - f_1 \zeta _{tt} - g\right ), \end{equation}

where $p$ denotes the pressure and $\rho$ the fluid density.

The equation is vertically integrated from an arbitrary elevation, $z$ , to the free surface, retaining only the leading-order terms:

(4.6) \begin{equation} p\mid _{z}^{\zeta } = p_a - p(z) = \rho \left (\left (\int _z^0 f_{1_x} \textrm{d}z\right ) u_{0_t} - \left (\int _z^0 f_1 \textrm{d}z\right ) \zeta _{tt} - g(\zeta -z)\right ), \end{equation}

after which the pressure is isolated:

(4.7) \begin{equation} p(z) = p_a + \rho \left (g(\zeta -z) - \left (\int _z^0 f_{1_x} \textrm{d}z\right ) u_{0_t} + \left (\int _z^0 f_1 \textrm{d}z\right ) \zeta _{tt}\right ). \end{equation}

Taking the derivative in the $x$ direction:

(4.8) \begin{align} p_x &= \rho \left (g\zeta _x - \left (\int _z^0 f_{1_x} \textrm{d}z\right )_x u_{0_t} - \left (\int _z^0 f_{1_x} \textrm{d}z\right ) u_{0_{xt}} + \left (\int _z^0 f_1 \textrm{d}z\right )_x \zeta _{tt} \right.\notag \\ &\quad+ \left. \left (\int _z^0 f_1 \textrm{d}z\right ) \zeta _{ttx}\right ). \end{align}

Defining $e = \int _z^0 f \,\textrm{d}z$ , the Leibniz integration rule is utilised to derive

(4.9) \begin{equation} p_x = \rho \left (g\zeta _x - e_{xx} u_{0_t} -e_x u_{0_{xt}} +e_x\zeta _{tt} + e\zeta _{ttx}\right ). \end{equation}

This expression is substituted into the horizontal component of the momentum equation, applied in a vertically averaged sense:

(4.10) \begin{align} &\qquad\quad\int _{-h}^0 u_t \,\textrm{d}z = -\frac {1}{\rho }\int _{-h}^0 p_x\, \textrm{d}z ,\end{align}

(4.11) \begin{align} &\left (1-E_2)\right ) u_{0_t} = -g h\zeta _x - 2 E_1 \zeta _{tt} - E_0 \zeta _{ttx}, \end{align}

where $E_0 = \int _{-h}^0 e \, \textrm{d}z$ , $E_1 = \int _{-h}^0 e_{x}\, \textrm{d}z$ and $E_2 = \int _{-h}^0 e_{xx} \,\textrm{d}z$ .

We then divide by the left-hand coefficient and differentiate with respect to $x$ :

(4.12) \begin{align} u_{0_{xt}} = -\frac {\partial }{\partial x} \left (\frac {g h\zeta _x + 2 E_1 \zeta _{tt} + E_0 \zeta _{ttx}}{1-E_2}\right ). \end{align}

The linear relation $u_{0_{xt}} = -\zeta _{tt}$ is utilised to eliminate $u_0$ and provide an equation as a function of $\zeta$ :

(4.13) \begin{align} \zeta _{tt} = \frac {\partial }{\partial x} \left (\frac {g h\zeta _x + 2 E_1 \zeta _{tt} + E_0 \zeta _{ttx}}{1-E_2}\right ). \end{align}

This equation is referred to as the LTDPCE in the following sections.

4.2. The 3-D equation

The same procedure is applied here to the 3-D case. As in the equation in the frequency domain, a local axis system is oriented such that the $x$ axis is aligned with the bottom gradient and perpendicular to the $y$ axis. The same profiles are utilised:

(4.14) \begin{align} \mathbf{u}(x,y,z,t) &= \left (u(x,y,z,t),v(x,y,z,t) \right ) \notag\\&= \left (f_{1_z}(h,h',z) u_0(x,y,t), f_{2_z}(h,h',z) v_0(x,y,t) \right ), \end{align}

(4.15) \begin{align} &\qquad\boldsymbol{\Psi }(x,y,z,t) = \left (f_1 u_0,f_2 v_0 \right ), \end{align}

(4.16) \begin{align} w(x,&y,z,t) = -\nabla \cdot \boldsymbol{\Psi } = -f_{1_{x}} u_0 - f_1 u_{0_{x}} - f_2 v_{0_{y}}, \end{align}

where $f_1$ and $f_2$ are defined in equations (3.2) and (3.3).

The relation to the free-surface location is provided by the linear kinematic free-surface condition:

(4.17) \begin{equation} \zeta _t(x,y,z,t) = w\mid _{z=0} = -u_{0_{x}} - v_{0_{y}}. \end{equation}

The bottom boundary condition, $w + \nabla \cdot \mathbf{u} = 0$ , is satisfied in this representation.

Next, the horizontal component of the linearised momentum equation is utilised:

(4.18) \begin{equation} \rho w_t = -p_z - \rho g, \end{equation}

and by substituting the vertical velocity profile, we derive

(4.19) \begin{equation} p_z = \rho \left (f_{1_{x}} u_{0_t} + f_1 u_{0_{x t}} + f_2 v_{0_{y t}} - g\right ). \end{equation}

The equation is vertically integrated from an arbitrary elevation $z$ to the free surface while retaining only the leading-order terms:

(4.20) \begin{align} p\mid _{z}^{\zeta } &= p_a - p(x,y,z) = \rho \left (\left (\int _z^0 f_{1_{x}} \textrm{d}z\right ) u_{0_t} \right.\notag\\ &\quad +\left. \left (\int _z^0 f_1 \textrm{d}z\right ) u_{0_{x t}} + \left (\int _z^0 f_2 \textrm{d}z\right ) v_{0_{y t}} - g(\zeta -z)\right ), \end{align}

after which the pressure is isolated:

(4.21) \begin{equation} p = p_a + \rho \left (g(\zeta -z) -\left (\int _z^0 f_{1_{x}} \textrm{d}z\right ) u_{0_t} - \left (\int _z^0 f_1 \textrm{d}z\right ) u_{0_{x t}} - \left (\int _z^0 f_2 \textrm{d}z\right ) v_{0_{y t}} \right ). \end{equation}

Taking the derivatives in the horizontal directions:

(4.22) \begin{align} p_{x} &= \rho \left (g\zeta _{x} -\left (\int _z^0 f_{1_{xx}} \textrm{d}z\right ) u_{0_t} - 2 \left (\int _z^0 f_{1_{x}} \textrm{d}z\right ) u_{0_{x t}} - \left (\int _z^0 f_1 \textrm{d}z\right ) u_{0_{xxt}} \right.\notag\\&\quad - \left.\left (\int _z^0 f_{2_{x}} \textrm{d}z\right ) v_{0_{yt}} - \left (\int _z^0 f_2 \textrm{d}z\right ) v_{0_{xyt}} \right ), \end{align}

(4.23) \begin{align} p_{y} &= \rho \left (g\zeta _{y} - \left (\int _z^0 f_{1_{x}} \textrm{d}z\right ) u_{0_{yt}} - \left (\int _z^0 f_1 \textrm{d}z\right ) u_{0_{xyt}} - \left (\int _z^0 f_2 \textrm{d}z\right ) v_{0_{yyt}} \right ). \end{align}

Defining $g_1 = \int _z^0 f_1 \,\textrm{d}z$ , $g_2 = \int _z^0 f_2\, \textrm{d}z$ and utilising the Leibniz integration rule:

(4.24) \begin{align} &\qquad\qquad p_{y} = \rho \left (g\zeta _{y} - g_{1_{x}} u_{0_{yt}} - g_{1} u_{0_{xyt}} - g_2 v_{0_{yyt}} \right ), \end{align}

(4.25) \begin{align} &p_{x} = \rho \left (g\zeta _{x} - g_{1_{xx}} u_{0_t} - 2 g_{1_{x}} u_{0_{xt}} - g_1 u_{0_{xxt}} - g_{2_{x}} v_{0_{yt}} - g_2 v_{0_{xyt}} \right ). \end{align}

These expressions are now substituted into the horizontal components of the momentum equation, applied in a vertically averaged sense:

(4.26) \begin{align} &F_1 u_{0_t} = - g h \zeta _{x} + G_{1c} u_{0_t} + 2 G_{1b} u_{0_{xt}} + G_{1a} u_{0_{xxt}} + G_{2b} v_{0_{yt}} + G_{2a} v_{0_{xyt}}, \\[-10pt] \nonumber \end{align}

(4.27) \begin{align} &\qquad\qquad F_2 v_{0_t} = - g h \zeta _{y} + G_{1b} u_{0_{yt}} + 2 G_{1a} u_{0_{xyt}} + G_{2a} v_{0_{yyt}}, \\[10pt] \nonumber \end{align}

where $F_1 = \int _{-h}^0 f_{1_z} \textrm{d}z$ , $F_2 = \int _{-h}^0 f_{2_z} \textrm{d}z$ , $G_{1a} = \int _{-h}^0 g_1 \textrm{d}z$ , $G_{1b} = \int _{-h}^0 g_{1_{x}} \textrm{d}z$ , $G_{1c} = \int _{-h}^0 g_{1_{xx}} \textrm{d}z$ , $G_{2a} = \int _{-h}^0 g_2 \textrm{d}z$ and $G_{2b} = \int _{-h}^0 g_{2_x} \textrm{d}z$ .

Equations (4.17), (4.26) and (4.27) provide a system of partial differential equations that can be solved for the vector stream function components and free-surface location. The equations are rotated to the fixed $x,y$ axis system in the same way as in the 3-D PCSSE.

5. Simulations

Numerical simulations of the LPCSSE and the LTDPCE are compared with those of the ESWE, as these equations provide a first-order long-wave approximation. Additionally, the PCSSE derived by Schwartz et al. (Reference Schwartz, Oron and Agnon2023) was simulated, with no visible differences observed compared with the LPCSSE in the examples and regions considered in the presented simulations. The simulations are achieved by numerically calculating the coefficients and using the NDsolve function in Mathematica. All simulations are time-harmonic, ensuring no additional treatment is needed for the time-dependent models. The spatial step size is approximately 0.01. They are conducted on several bottom topographies (figure 4), including Booij’s ramp, a constant-slope beach with normal and oblique incident waves, Roseau’s bathymetry and a class II Bragg resonance. The results are presented in the following figures.

Figure 4. Illustration of the test case bathymetries.

Figure 5. Reflection coefficient for a progressive wave flowing over Booij’s ramp. (a) The full domain. (b) Enlargement of a part of the region.

Figure 6. Normalised free-surface values versus $\sigma h$ for a normal incidence simulation. Regular standing wave with (a) $45^\circ$ bed slope and (b) $30^\circ$ bed slope.

The bed shape in the form of a ramp, as tested by Booij (Reference Booij1983), is simulated in figure 5. The ramp consists of a constant-slope region between two semi-infinite horizontal segments. Both PC models provide a better match to the exact solution in the range $\sigma h_0$ between $8$ and $15$ as shown in figure 5(b), as well as for waves shorter than $\sigma h_0 = 4$ . For longer waves, the differences become indistinguishable.

The next test case involves normal and oblique incident waves on a constant-slope beach. An exact solution for waves that are normally incident to the bottom gradient is found in Stoker (Reference Stoker1957) for bottom slopes of ${\varPi }/{2n}$ , where $n$ is an integer. A propagating wave solution is given as the sum of two standing waves: a regular standing wave, finite at the origin, and a singular standing wave, infinite at the origin. Simulations of the regular standing wave were performed on bed slopes up to $45^{\circ}$ . A relative error $\displaystyle E_r = ({R_m-R_a})/({R_a})$ is defined, where $R_m$ and $R_a$ are the PCMSE/ESWE model and the analytical reflection coefficients, respectively. For $45^{\circ}$ degrees (figure 6 a) the LPCSSE provides a very accurate result up to $\sigma h_0 = 1.3$ ( $E_r = 0.0015$ ), while the ESWE deviates from the exact solution well before that ( $E_r = 0.17$ at $\sigma h_0 = 1.3$ ). The LTDPCE closely follows the LPCSSE until $\sigma h_0 = 1.2$ , then remains between the LPCMSE and ESWE ( $E_r = 0.025$ at $\sigma h_0 = 1.3$ ). For $30^{\circ}$ degrees (figure 6 b), the results are similar but less pronounced. The analytical LPCSSE was simulated in this case, requiring only $App_1$ (2.8) for the calculation, and produced results identical to the approximated analytical equation.

Ehrenmark (Reference Ehrenmark1998) presented an exact solution for obliquely incident waves over a plane beach. The solution is limited, as in the normal incidence of Stoker, to bottom slopes of ${\varPi }/{2n}$ , where $n$ is an integer. Figure 7 shows the case of $45^{\circ}$ oblique incidence on a $45^{\circ}$ bottom slope. The ESWE model performs better up to $\sigma h = 1.2$ (maximum error of $0.034$ for LPCMSE, $0.056$ for LTDPCE and $0.007$ for ESWE), but beyond that, the PC equations provide more accurate results (for example, at $\sigma h = 1.5$ , the relative error is $0.064$ for LPCMSE, $0.011$ for LTDPCE and $0.18$ for ESWE).

Figure 7. Normalised free-surface values versus $\sigma h$ for an oblique incidence simulation and an incidence angle of $45^\circ$ .

Roseau (Reference Roseau1976) derived the only analytical solution for waves propagating over a non-constant-slope bathymetry. In his solution, the bed location is defined by

(5.1) \begin{align} &\frac {x(\xi )}{h_0} = \xi - (2 \varPi \beta )^{-1} \left (1-\frac {h_L}{h_0}\right ) \ln \left (1+\textrm{e}^{2\beta \varPi \xi } + 2\textrm{e}^{\beta \varPi \xi } \cos(\beta \varPi ) \right ), \end{align}

(5.2) \begin{align} &\quad\frac {z(\xi )}{h_0} = -1 + (\varPi \beta )^{-1} \left (1-\frac {h_L}{h_0}\right ) \arctan \left (\frac {\sin (\beta \varPi )}{\textrm{e}^{-\beta \varPi \xi } +\cos (\beta \varPi )} \right ), \end{align}

where $\beta \in (0,1)$ is a shoaling parameter and the bed is defined parametrically as $z(\xi ) = -h (x(\xi ) )$ . The bed function tends to two flat semi-infinite sections at the limits $x \to \pm \infty$ . In figure 8, a comparison of the reflection coefficients for the steepness value $\beta = 0.5$ is presented. The LPCSSE provides a significantly better match than the ESWE up to approximately $\sigma h_0 = 1.2$ . The LTDPCE is even better, with results very close to the exact solution up to $\sigma h_0 = 1.6$ . This is the region where the long-wave equations are expected to perform. The relative error of the reflection coefficients is presented in figure 8(b). For more moderate slopes, i.e. for smaller $\beta$ values, the errors of all models reduce significantly.

Figure 8. (a) Simulation of the reflection coefficient versus $\sigma h_0$ for the Roseau bathymetry; $\beta = 0.5$ , $ {h_L}/{h_0} = 0.1$ . (b) Simulation of the relative error versus $\sigma h_0$ .

Figure 9. Comparison of the reflection coefficient values obtained from simulations of the LPCSSE (red dot-dashed curve), ESWE (blue small-dashed curve), LPCTDE (orange large-dashed curve) and the experimental results obtained by Guazzelli, Rey & Belzons (Reference Guazzelli, Rey and Belzons1992) as shown in figure 2 of that article (black solid curve).

The classical MSE, which accounts for second-order terms, generally provides a good agreement with the class I Bragg resonance. This phenomenon involves the interaction of two surface waves and one bottom undulation that satisfy predefined conditions (see e.g. Liu & Yue Reference Liu and Yue1998). This is not the case in class II, for which they fail to provide a reasonable approximation. Class II involves the interaction of two surface waves and two bottom undulations. A simulation of the impact of class II on the free-surface values is provided in figure 9. This test case was presented by Guazelli et al. (Reference Guazzelli, Rey and Belzons1992) and involves two bed wavenumbers, $k_{b_1} = {\pi }/{6} \, \textrm{cm}^{-1}$ and $k_{b_2} = {\pi }/{3} \, \textrm{cm}^{-1}$ , a patch length of $L = 48 \, \textrm{cm}$ outside of which the bed is horizontal and a ratio of the average bottom depth to the bed wavenumber amplitude of ${\Delta H}/{H_0} = 0.25$ and $H_0 = 2.5\, \textrm{cm}$ . Guazelli’s model, which is used for the comparison, accounts for three evanescent modes. In the first peak ( $\sigma h \approx 1.7$ ), the PC-based models show a significantly better match than the ESWE compared with the numerical results of Guazelli. This peak is the result of the class II resonance. In the second peak ( $\sigma h \approx 3.1$ ), the LTDPCE and ESWE are very far off. The LPCSSE shows a reasonable match with the peak location and an excellent match with the width and maximum value of the peak curve. This peak is the result of class I resonance.