2.1. The angular profile

2.1.1. Cylindrical coordinates – problem definition

The governing equations and boundary conditions are presented in cylindrical coordinates in order to provide an approximate angular profile for the two-component vector streamfunction in a 3-D fluid domain. This domain is contained in the polar plane for which $r$ and $\theta$ are the radial distance and the angle, respectively, and the longitudinal $y$-axis is normal to the polar plane. The ray $\theta =0$ represents the mean location of the free surface, and $r=0$ corresponds to the shore, as depicted in figure 1.

Figure 1. Polar coordinates – problem description.

Introducing the complex-valued vector streamfunction

(2.1)\begin{equation} \boldsymbol{\boldsymbol{\varPsi}}(x,y,z) \equiv \int_{{-}h}^z \boldsymbol{u}(x,y,z_0) \,{\rm d} z_0, \end{equation}

which defines the time harmonic streamfunction by ${\rm Re}(\boldsymbol {\varPsi } \exp ({-{\rm i}\omega t}))$, with $\omega$ being the angular frequency of the incoming wave, and $z=-h$ being the location of the bottom. The horizontal and vertical velocity components are denoted by $\boldsymbol {u}$ and $w$, respectively, and are determined via

(2.2a,b)\begin{equation} \boldsymbol{u} = \frac{\partial \boldsymbol{\varPsi}}{\partial z},\quad w ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varPsi}, \end{equation}

with $\boldsymbol {\nabla } = ({\partial }/{\partial x}, {\partial }/{\partial y})$. The vector streamfunction $\boldsymbol {\varPsi }$ satisfies the Laplace equation formulated in cylindrical coordinates in the flow domain

(2.3)\begin{gather} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial \boldsymbol{\varPsi}}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 \boldsymbol{\varPsi}}{\partial \theta^2}+ \frac{\partial^2 \boldsymbol{\varPsi}}{\partial y^2}= 0 \quad {\rm in}\ 0< r<\infty,\quad -\alpha<\theta<0,\enspace -\infty< y<\infty. \end{gather}

The bottom boundary condition is formulated in terms of the vector streamfunction (2.1)

(2.4)\begin{equation} \boldsymbol{\varPsi} = 0 \quad {\rm at}\ \theta ={-}\alpha, \end{equation}

and the Cartesian free-surface condition is translated to cylindrical coordinates. At the undisturbed free surface

(2.5a,b)\begin{equation} r {\rm d} \theta = {\rm d} z \rightarrow \frac{1}{r} \frac{\partial \boldsymbol{\varPsi}}{\partial \theta} = \frac{\partial \boldsymbol{\varPsi}}{\partial z},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varPsi} = \boldsymbol{\nabla}_{cy} \boldsymbol{\cdot} \boldsymbol{\varPsi}, \end{equation}

with $\boldsymbol {\nabla }_{cy} = ({\partial }/{\partial r},{\partial }/{\partial y})$.

These relations are then substituted into the boundary condition to be satisfied at the undisturbed location $z=0$ of the free surface. These conditions arise from the minimization of the energy functional (Kim & Bai Reference Kim and Bai2004) in terms of the vector streamfunction $\boldsymbol {\varPsi }$. They are then written in the Cartesian frame in order to derive the equivalent cylindrical condition

(2.6)\begin{equation} \frac{\partial \boldsymbol{\varPsi}}{\partial z} +\frac{1}{\sigma} \boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\varPsi}) = 0 \quad {\rm at}\ z = 0, \end{equation}

where $\sigma \equiv {\omega ^2}/{g}$. In cylindrical coordinates, (2.6) is written as

(2.7)\begin{equation} \frac{1}{r} \frac{\partial \boldsymbol{\varPsi}}{\partial \theta} + \frac{1}{\sigma} \boldsymbol{\nabla}_{cy} (\boldsymbol{\nabla}_{cy} \boldsymbol{\cdot} \boldsymbol{\varPsi})= 0 \quad {\rm at} \theta = 0. \end{equation}

Equations (2.3)–(2.7) constitute a closed set of equations to find a solution for the vector streamfunction $\boldsymbol {\varPsi }$, when supplemented by lateral boundary conditions. This solution can then be utilized to determine the free-surface location by resolving it from the kinematic boundary condition

(2.8)\begin{equation} \eta = \frac{1}{{\rm i} \omega} \boldsymbol{\nabla}_{cy} \boldsymbol{\cdot} \boldsymbol{\varPsi} \quad {\rm at}\ \theta = 0. \end{equation}

2.1.2. The $\theta$-dependent profile

The set of equations (2.3)–(2.7) is solved by assuming a solution to the vector streamfunction, $\boldsymbol {\varPsi }$, in a separable form. It is noted that this variable separation will provide an approximate solution to the problem when solved on the cylinder segment. The segment is composed of a circle sector in the polar plane, defined between the free surface and the planar bottom and along the shore in the longitudinal plane

(2.9)\begin{equation} \boldsymbol{\varPsi}(r,\theta,y) = \sum_{n=0}^{\infty} R_n(r) F_n(\theta) Y_n(y), \end{equation}

and $F_n(\theta )$ will be utilized as the approximation for the angular profile of the vector streamfunction.

Substituting this expression into (2.3) yields

(2.10)\begin{equation} \left.\begin{gathered} \frac{1}{F_n} \frac{{\rm d}^2 F_n}{{\rm d} \theta^2} = \gamma_n^2,\\ \frac{1}{Y_n} \frac{{\rm d}^2 Y_n}{{{\rm d}y}^2} ={-}\delta_n^2,\\ \frac{{\rm d}^2 R_n}{{\rm d} r^2} + \frac{1}{r} \frac{{\rm d} R_n}{{\rm d} r} + \left(\frac{\gamma_n^2}{r^2} - \delta_n^2\right)R_n = 0, \end{gathered}\right\} \end{equation}

where $\gamma _n$ and $\delta _n$ are positive constants, analogous to the components of the vector wavenumber $k$, normally used in the Cartesian formulation.

Defining the solution around a local wavenumber expressed by the constants $\gamma _n$ and $\delta _n$, an oscillating solution is chosen in the $y$ direction. In the $r$ direction, Bessel functions are obtained. The solution in the angular direction has exponential form

(2.11)\begin{equation} F(\theta) = A \sinh({\gamma \theta}) + B \cosh({\gamma \theta}). \end{equation}

This solution is then substituted into the boundary condition at the bottom and then normalized with the value at the free surface ($\theta =0$) to yield

(2.12)\begin{equation} F(\theta) = \frac{\sinh(\gamma(\theta+\alpha))}{\sinh(\gamma \alpha)}. \end{equation}

The $\sinh (\gamma (\theta +\alpha ))$ shape will serve in what follows as an approximated profile and will be utilized to derive the mild-slope equation in the horizontal direction.

This profile is then substituted into the time-averaged cylindrical Lagrangian in order to provide a solution for the case of a 3-D planar beach

(2.13)\begin{equation} \left.\begin{gathered} \bar{L} = \int \int L \,{\rm d} r \,{{\rm d}y},\\ \bar{L} = \frac{\rho}{2} \int_{-\alpha}^0 \left(|\boldsymbol{\nabla}_{cy} \boldsymbol{\cdot} \boldsymbol{\varPsi}|^2 + \frac{1}{r^2} \left|\frac{\partial \boldsymbol{\varPsi}}{\partial \theta}\right|^2 \right) r \,{\rm d}\theta - \frac{\rho}{2 \sigma} |\boldsymbol{\nabla}_{cy} \boldsymbol{\cdot} \boldsymbol{\varPsi}(r,y,0)|^2. \end{gathered}\right\} \end{equation}

Here, $L$ is the Lagrangian density and the double integration is made in the domain defined by lateral boundary conditions which appropriately depict the problem at hand. Since the wave field is harmonic, the (quadratic) Lagrangian has a zeroth harmonic component and a second harmonic one. Hence, there are two alternative applications of the Lagrangian: the present, time-averaged component, utilized by Kim & Bai (Reference Kim and Bai2004), or the oscillatory component, utilized by Chamberlain & Porter (Reference Chamberlain and Porter1995).

The extension of this formulation to a full 3-D bathymetry will be provided in future work.

For the case of a general 2-D bathymetry, the Cartesian equivalent to the angular profile is substituted into the Cartesian Lagrangian

(2.14)\begin{equation} \left.\begin{gathered} \bar{L} = \int L {{\rm d}x},\\ L = \frac{\rho}{2} \int_{{-}h(x)}^0 \left(\left|\frac{\partial \boldsymbol{\varPsi}}{\partial x}\right|^2 + \left|\frac{\partial \boldsymbol{\varPsi}}{\partial z}\right|^2 \right) {\rm d} z - \frac{\rho}{2 \sigma} \left|\frac{\partial \boldsymbol{\varPsi}(x,0)}{\partial x}\right|^2, \end{gathered}\right\} \end{equation}

after which the first variation is equated to zero to derive a 2-D MSE.

2.2. The 2-D PC mild-slope equation

In two dimensions, the problem is constructed in the polar plane and the vector streamfunction reduces to a scalar streamfunction $\boldsymbol {\varPsi } = (\psi,0)$. The set of equations (2.3)–(2.7) reduces to

(2.15)\begin{gather} \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 {\rm in}\ 0< r<\infty,\quad -\alpha<\theta<0. \end{gather}

(2.16)\begin{gather}\psi = 0 \quad {\rm at}\ \theta ={-}\alpha, \end{gather}

(2.17)\begin{gather}\frac{1}{r} \frac{\partial \psi}{\partial \theta} +\frac{1}{\sigma} \frac{\partial^2 \psi}{\partial r^2}= 0 \quad {\rm at}\ \theta = 0. \end{gather}

The solution form for variable separation around a local wavenumber

(2.18)\begin{gather} F(\theta) = A \sinh(\gamma {\theta}) + B \cosh({\gamma \theta}), \end{gather}

(2.19)\begin{gather}R(r) = C \cos(\gamma \ln r) + D \sin(\gamma \ln r), \end{gather}

where the constants are determined by substitution into the boundary conditions. A localized dispersion relation is derived when substituting the solution form into the free-surface boundary condition (2.17)

(2.20)\begin{equation} \sigma r = \gamma \tanh(\gamma\alpha). \end{equation}

This is simply a reconstruction of the relation derived by Ehrenmark (Reference Ehrenmark2005), as stated above.

When solving for a general bathymetry, the local bed gradient defines the angle $\alpha$. The origin of the polar reference frame is located where the corresponding ray intersects with the line of the undisturbed free surface. Therefore, at each local point the axis origin shifts according to the bottom gradient.

The relations between the polar and Cartesian parameters, as depicted in figure 2, are provided by the following:

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

The following form of $\psi$ is substituted into the Lagrangian

(2.22)\begin{equation} \psi(r,\theta) = \psi_0(r) F(\theta), \end{equation}

and the problem is then to be solved for the unknown function $\psi _0(r)$.

Figure 2. Polar-Cartesian relations.

The Cartesian equivalent profile of (2.12) is given by

(2.23)\begin{equation} F \left(h(x),h'(x),z \right) = \frac{\sinh\left(\gamma \left(\arctan\left(\dfrac{z h'}{h}\right) + \arctan(h') \right)\right)}{\sinh \left(\gamma \arctan(h') \right)}. \end{equation}

A solution for $\psi (x,z)$ is sought in the form of

(2.24)\begin{equation} \psi(x,z) = \psi_0(x) F \left(h(x),h'(x),z \right). \end{equation}

Substituting the profile into the Lagrangian (2.14), yields

(2.25)\begin{equation} L = \left(a - \frac{1}{\sigma} \right) | \psi_{0_x} |^2 + (c+d) | \psi_0 |^2 + 2b \,{\rm Re} \left(\psi_0 \psi_{0_x}^* \right), \end{equation}

where $\textrm {Re}(x)$ denotes the real part of $x$, $*$ denotes the complex conjugate and the coefficients $a,b,c$ and $d$ are real-valued functions of $h(x)$ and $h'(x)$ defined by

(2.26a–d)\begin{equation} a = \int_{{-}h}^0 F^2 \,{\rm d} z \quad b = \int_{{-}h}^0 F F_x \,{\rm d} z \quad c = \int_{{-}h}^0 F_x^2 \,{\rm d} z \quad d = \int_{{-}h}^0 F_z^2 \,{\rm d} z. \end{equation}

Following the substitution, the Lagrangian becomes a functional of a single variable, $\psi _0$, and the Euler–Lagrange equation is derived in the form

(2.27)\begin{equation} -\left[ \left(a-\frac{1}{\sigma}\right) \psi_{0_{x}} + b \psi_{0} \right]_x + b \psi_{0_x} + (c+d) \psi_0 = 0. \end{equation}

Its solution will provide the stationary point of the Lagrangian and the solution to the original system of equations. Note that the coefficients of (2.27) depend on $x$, so the $b$-term does not cancel out.

This is the PCMSE and the main result presented in the paper. The structure of the equation is similar to that of the CMSE and the fundamental difference is in the order of the bottom derivative in the coefficients. This is up to third order here as opposed to second order in the CMSE.

We note that numerical issues arise with evaluation of the coefficients $a,b$ and $c$ defined in (2.26a–d) when solving (2.27) in a horizontal or close to horizontal bed. The values of $\alpha$ and $\theta$ both tend to zero at the limit where $h^\prime \to 0$. In order to obtain a smooth solution in the entire range from a steep slope to a horizontal bottom, the variables are normalized in the following form:

(2.28a–d)\begin{equation} \alpha_0 = \alpha r_0,\quad \theta_0 = \theta r_0,\quad \gamma_0 = \frac{\gamma}{r_0},\quad r_0 \equiv \frac{h}{h^\prime}. \end{equation}

Under this transformation, the vertical velocity profile remains unchanged. The polar coordinate variables reduce to the equivalent Cartesian variables in the limit of a horizontal bottom, and this resolves the numerical complications

(2.29a,b)\begin{equation} \lim_{h^\prime \to 0} \alpha_0 = h,\quad \lim_{h^\prime \to 0} \theta_0 = z, \end{equation}

and the polar dispersion relation reduces to its Cartesian version

(2.30a–c)\begin{equation} \sigma = \gamma_0 \tanh(\mu \gamma_0 h),\quad \mu \equiv \alpha \cot(\alpha),\quad \lim_{h^\prime \to 0} \mu = 1. \end{equation}

This representation reduces the PCMSE to the CMSE for bottom slopes which tend to zero. Effectively, this means that when solving horizontal or close to horizontal bottom profiles the CMSE and PCMSE will provide similar results. For relatively steep gradients the results will differ according to the coefficient terms of each equation.

In terms of the normalized variables, the streamfunction profile is now given as

(2.31)\begin{equation} F(\gamma_0,h,h^\prime) = \frac{\sinh\left[\gamma_0(\theta_0+\alpha_0)\right]}{\sinh(\gamma_0\alpha_0)} = \frac{\sinh\left[\dfrac{\gamma_0 h}{h^\prime}\left(\arctan\left(\dfrac{z h^\prime}{h}\right) + \arctan(h^\prime) \right)\right]}{\sinh \left(\dfrac{\gamma_0 h}{h^\prime} \arctan(h^\prime) \right)}. \end{equation}

This profile and its derivatives are integrated with respect to the vertical coordinate in order to obtain the coefficients of (2.27). Due to the complexity of the terms derived by analytical integration, this task is accomplished in two ways. The first one is by a Taylor expansion, assuming a small $\alpha$, which reduces the complexity of the terms derived from the analytical integration. The challenge is that using $\alpha$ as the expansion parameter requires taking several orders to obtain values that warrant a satisfactory approximation. This leads to lengthy and cumbersome expressions. To overcome this difficulty and to simplify the coefficients resulting from (2.26a–d), the following parameter is utilized:

(2.32)\begin{equation} p \equiv 1-\mu. \end{equation}

With $p$ as the expansion parameter, accurate simulations are derived by expanding (2.31) to first order. Also, $\alpha$ is expanded into powers of $p$ by using successive approximations $\alpha = \sqrt {3 p} (1-{p}/{10} + O(p^2) )$ and subsequently used in the integration (2.26a–d) to obtain the coefficients of the PCMSE. The coefficient terms to first order in $p$ are provided in Appendix A. This is referred to in the following as the PCMSE.

The second is by numerical integration of the full coefficients in what will be referred to as the ‘steep-slope’ equation. In this equation, no approximations have been imposed on the bottom slope and therefore it exceeds the mild-slope approximation and provides a first ‘steep-slope’ equation. This is referred to in the following as the PCSSE (PC steep-slope equation).

2.3. The quasi-3-D equation of oblique incidence

The case of an incoming wave which is incident to the bottom gradient at a varying angle is considered. A local axis system is defined such that the $x$ axis is oriented in the direction of the gradient at each point. The PC vertical profile is utilized in this direction of varying depth topography (2.31). This profile is denoted as $f_1$ in the following derivation. The corresponding local $y$ axis is that of a horizontal bottom which specifies a harmonic solution dependence. The profile in this direction is defined as

(2.33)\begin{equation} f_2 = \frac{\sinh\left(\gamma_0(z+h)\right)}{\sinh(\gamma_0 h)}. \end{equation}

Where $\gamma _0$ is the solution of the PC dispersion relation (2.30a–c), as stated previously. The vector streamfunction is approximated by a separation of variables. This is expressed as a multiplication of a $y$ dependent harmonic function, an approximated vertical profile and an unknown horizontal function. The horizontal function represents the stream value at the undisturbed free surface

(2.34)\begin{equation} \boldsymbol{\varPsi}(x,y,z) = \exp({{\rm i} \beta_0 y}) \left(f_1(h,h',z) \psi_1(x), f_2(h,z) \psi_2(x)\right), \end{equation}

with $\beta _0$ being the $y$ component of the wavenumber vector. This is then substituted into the 3-D time-averaged Lagrangian

(2.35)\begin{equation} L = \int_{{-}h}^0 \left( | u |^2 + | v |^2 + | w |^2 \right) {\rm d} z - g |\zeta|^2, \end{equation}

where $u,v$ and $w$ are the complex-valued local velocities in the $x,y$ and $z$ directions and $\zeta$ is the complex-valued free-surface location. The velocities are calculated according to the streamfunction definition

(2.36)\begin{equation} \left.\begin{gathered} {}[u,v] = \boldsymbol{\varPsi}_z = \exp({{\rm i} \beta_0 y}) \left(f_{1_z} \psi_1, f_{2_z} \psi_2 \right)\\ w ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varPsi} ={-}\exp({{\rm i} \beta_0 y}) (f_{1_x} \psi_1 + f_1 \psi_{1_x} + {\rm i} \beta_0 f_2 \psi_2 ), \end{gathered}\right\} \end{equation}

and the free-surface location is calculated according to the kinematic linear free-surface condition

(2.37)\begin{equation} \zeta = \frac{1}{{\rm i}\omega} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\varPsi}_0 = \frac{1}{{\rm i} \omega} \exp({{\rm i} \beta_0 y})(\psi_{1_x} + {\rm i}\beta_0 \psi_2). \end{equation}

These expressions are substituted into the Lagrangian

(2.38)\begin{align} L &= (c_1+d_1) | \psi_1 |^2 + a_2 (\beta_0^2 - \gamma_0^2) | \psi_2 |^2 + a_1 | \psi_{1_x} |^2 \nonumber\\ &\quad + 2 b_1 (\psi_1 \psi_{1_x}^* + \psi_1^* \psi_{1_x}) + {\rm i} \beta_0 b_{12} (\psi_2 \psi_1^* - \psi_1\psi_2^*) \nonumber\\ &\quad + {\rm i}\beta_0 a_{12} (\psi_2 \psi_{1_x}^* - \psi_{1_x} \psi_2^*), \end{align}

where

(2.39)\begin{equation} \left.\begin{gathered} a_1 = \int_{{-}h}^0 f_1^2 \,{\rm d} z - \frac{1}{\sigma} \quad b_1 = \int_{{-}h}^0 f_1 f_{1_x} \,{\rm d} z \quad c_1 = \int_{{-}h}^0 f_{1_x}^2 \,{\rm d} z\\ d_1 = \int_{{-}h}^0 f_{1_z}^2 \,{\rm d} z \quad a_2 = \int_{{-}h}^0 f_2^2 \,{\rm d} z - \frac{1}{\sigma} \quad a_{12} = \int_{{-}h}^0 f_1 f_2 \,{\rm d} z - \frac{1}{\sigma}\\ b_{12} = \int_{{-}h}^0 f_2 f_{1_x} \,{\rm d} z. \end{gathered}\right\} \end{equation}

Variations are taken with respect to both parameters, $\psi _1$ and $\psi _2$, and equated to 0. The variation with respect to $\psi _2$ provides

(2.40)\begin{equation} \psi_2 = \frac{{\rm i} \beta_0 \left( b_{12} \psi_1 + a_{12} \psi_{1_x} \right)}{a_2 (\beta_0^2 - \gamma_0^2)} \equiv {\rm i}\beta_0 (r_1 \psi_1+r_2 \psi_{1_x}), \end{equation}

and with respect to $\psi _1$, after substituting for $\psi _2$ and rearranging

(2.41)\begin{equation} \psi_{1_{xx}} + \frac{R_1}{L_1} \psi_{1_x} + \frac{R_2}{L_1} \psi_1 = 0 \end{equation}

where

(2.42)\begin{equation} \left.\begin{gathered} R_1 ={-}a_{1_x}-\beta_0^2(b_{12}-a_{{12}_x}) r_2 + \beta_0^2 a_{12} (r_1+r_{2_x})\\ R_2 = c_1+d_1-b_{1_x}-\beta_0^2(b_{12}-a_{{12}_x}) r_2 + \beta_0^2 a_{12} (r_1+r_{2_x})\\ L_1 ={-}a_1+\beta_0^2 a_{12} r_2. \end{gathered}\right\} \end{equation}

3.1. Roseau's bathymetry

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

(3.1)\begin{gather} \left.\begin{gathered} \frac{x(\xi)}{h_0} = \xi - (2 \varPi \beta)^{{-}1} \left(1-\frac{h_L}{h_0}\right) \ln \left(1+\exp({2\beta \varPi \xi}) + 2\exp({\beta \varPi \xi}) \cos(\beta \varPi) \right),\\ \frac{z(\xi)}{h_0} ={-}1 + (\varPi \beta)^{{-}1} \left(1-\frac{h_L}{h_0}\right) \arctan \left(\frac{\sin(\beta \varPi)}{\exp({-\beta \varPi \xi}) +\cos(\beta \varPi)} \right), \end{gathered}\right\} \end{gather}

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$, as depicted in figure 3.

Figure 3. Roseau's bed function for ${h_L}/{h_0}=0.25$. The $x$ axis is the horizontal distance. The $z$ axis is the vertical location of the bottom. The values are normalized with respect to $h_0$. Green (dotted) curve – $\beta = 0.5$; red (dashed) curve – $\beta = 0.45$; blue (solid) curve – $\beta = 0.4$.

The reflection coefficient, defined as the ratio of the amplitude of a reflected wave to that of a corresponding incident wave, is as given by Roseau

(3.2)\begin{equation} | R| = \left| \frac{\sinh\left((k_0h_0-k_Lh_L)/\beta \right)}{\sinh\left((k_0h_0+k_Lh_L)/\beta \right)} \right|, \end{equation}

where $k_0$ and $k_L$ are the wavenumbers at the constant depth locations, $h_0$ and $h_L$, respectively. A solution to the case of reflection and transmission of temporally harmonic waves as presented by Porter (Reference Porter2019) is utilized here. The waves are of frequency $\omega$, incident from $x=-\infty$ over a finite region of variable 2-D bathymetry between two flat semi-infinite sections defined by $h = h_0\ \textrm {at}\ x<0$ and $h = h_L\ \textrm {at}\ x>L$,

This method is applied to the MSE models and the reflection coefficient is solved for numerically. The results are presented in comparison with the exact solution of Roseau for three values of $\beta$ and two relations between the depths at the flat regions ${h_L}/{h_0}$. The maximum slopes of the simulations which are presented here are $0.75,0.68$ and $0.62$ for $\beta = 0.5,0.45$ and $0.4$, respectively.

In addition, the conservation of energy flux (as presented by Massel Reference Massel1993) is tested for the steepest bathymetry case by equating net energy flux values on both sides of the variable region. This is stated by the following relation between the reflection and transmission coefficients and the relevant group velocity:

(3.3)\begin{equation} (C_g)_{h = h_0} (1-| R |^2) = (C_g)_{h = h_L} | T |^2 .\end{equation}

Figure 4 displays the reflection coefficient $R$ as a function of the non-dimensional $\sigma h_0$ for $\beta = 0.5$ for the CMSE, MMSE, PCMSE and the PCSSE in comparison with the analytical result of Roseau. In panel (b), $({h_l}/{h_0} = 0.25$) the results obtained from the PC equations are almost identical to the exact result and outperform the other models. In panel (a), (${h_l}/{h_0} = 0.1$) the differences are more significant and clearly the PC equations provide a better match than the other model equations. The PCSSE outperforms the PCMSE by a significant amount and is very close to the exact solution. In figure 5, the relative error ($E_r = {(R_m-R_r)}/{R_r}$, where $R_m$ and $R_r$ are the MSE model and Roseau reflection coefficients, respectively) is presented for the PCMSE and PCSSE in comparison with the MMSE and CMSE (a,c) and with the MSE with the standard Cartesian dispersion relation and with Ehrenmark's dispersion relation (b,d) for $\beta = 0.5$ as well. The energy flux conservation (3.3) was tested for the case of ${h_L}/{h_0} = 0.25$. The maximum relative income–outcome energy flux value

(3.4)\begin{equation} \left(E_f = \frac{(C_g)_{h = h_0} (1-| R |^2) - (C_g)_{h = h_L} | T |^2}{(C_g)_{h = h_0} (1-| R |^2) + (C_g)_{h = h_L} | T |^2}\right), \end{equation}

is $4\times 10^{-7}$ for all $\sigma h_0$ values which were tested. Figure 6 compares the relative error of the CMSE, MMSE and the PCMSE for $\beta$ values of $0.45$ and $0.4$. In panel (c), the PCMSE and the CMSE models provide indistinguishable results. The PCMSE can be seen to provide smaller relative errors in comparison with all models in an overall observation.

Figure 4. Simulation of the reflection coefficient vs $\sigma h_0$ for the Roseau bathymetry. Here, $\beta = 0.5$. Panels show (a) ${h_L}/{h_0} = 0.1$, (b) ${h_L}/{h_0} = 0.25$. Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; green (dotted) curve – MMSE; black (solid) curve - exact solution.

Figure 5. Simulation of the relative error vs $\sigma h_0$ for the Roseau bathymetry. Here, $\beta = 0.5$. Panels show (a,b) ${h_L}/{h_0} = 0.1$, (c,d) ${h_L}/{h_0} = 0.25$. (a,c) Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; green (dotted) curve – MMSE. (b,d) Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – MSE; green (dotted) curve – MSE with Ehrenmark's dispersion.

Figure 6. Simulation of the relative error vs $\sigma h_0$ for the Roseau bathymetry. Panels show (a,c) ${h_L}/{h_0} = 0.1$, (b,d) ${h_L}/{h_0} = 0.25$, for (a,b) $\beta = 0.45$, (b,d) $\beta = 0.4$. Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; green (dotted) curve – MMSE.

3.2. Constant-slope bathymetry

An exact solution to waves normally incident to the bottom gradient can be found in Stoker (Reference Stoker1957). The solution is limited to bottom slopes of ${{\rm \pi} }/{2n}$, where $n$ is an integer. Simulations were performed for the bed slopes of $45,30,22.5,18$ and $15$ degrees. The analytical solution (Stoker Reference Stoker1957) is compared with the results of numerical simulations for the PCMSE, CMSE and a full polar MSE in figure 7 for slopes of $45^{\circ }$ and $30^{\circ }$ degrees. For $45^{\circ }$, the polar equation agrees slightly better than the PCMSE with the analytical result, and both provide better agreement than the CMSE, as observed in the figure. For $30^{\circ }$, all models are effectively indistinguishable from the exact solution, as is for all angles smaller than $30^{\circ }$. This strengthens the result from the previous test case for which the PCMSE has an advantage for steep slopes. The polar model, which is specifically designed for this test case, has a small advantage over the PC equations. The results of the CMSE matched those of the PC equations for milder slopes.

Figure 7. Normalized free-surface values for a normal incidence simulation for two values of the bed slope. Panels show (a) $45^{\circ }$, (b) $30^{\circ }$. Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; green (dotted) curve – polar MSE; black (solid) curve – exact solution.

3.3. Booij's ramp

The bed shape in the form of a ramp tested by Booij (Reference Booij1983) is simulated in what follows. The bathymetry is given by

(3.5)\begin{equation} h(x) = h_m - \left(\frac{h_p-h_m}{L}\right) x,\quad 0< x< L, \end{equation}

where $h_m$ and $h_p$ are the constant depths upstream and downstream of the ramp, respectively. Following the studies by Kim & Bai (Reference Kim and Bai2004), we use the values

(3.6a,b)\begin{equation} \frac{h_p}{h_m} = \frac{1}{3},\quad \frac{\omega^2 h_m}{g} = 0.6, \end{equation}

with the non-dimensional ramp length ${\omega ^2 L}/{g}$ in the range between $0.1$ and $5$. In figure 8, the values of the reflection coefficient obtained using the PCMSE, CMSE and MMSE are compared with the result of the full linear theory of Bai & Yeung (Reference Bai and Yeung1974) using the localized finite-element method. The CMSE was solved in the domain containing the two points of a discontinuous bed slope by using the jump conditions suggested by Kim & Bai (Reference Kim and Bai2004). The PCMSE was solved using the same jump conditions extended to include the second-order slope discontinuity, as detailed below.

Figure 8. Comparison of the reflection coefficient derived using the exact solution and computed using various model equations for Booij's ramp. Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; green (dotted) curve – MMSE; black (solid) curve – exact solution.

The streamfunction solution over the horizontal bed on both sides of the ramp is provided by (43) and (44) in Kim & Bai (Reference Kim and Bai2004) as

(3.7)\begin{equation} \left.\begin{gathered} \psi_0(x) \equiv \psi_0^-(x) = \frac{\omega A}{k_m}(\exp({{\rm i}k_mx})-R\exp({-{\rm i}k_mx})),\quad x<0,\\ \psi_0(x) \equiv \psi_0^+(x) = \frac{\omega A}{k_p}T\exp({{\rm i}k_px}),\quad x>L, \end{gathered}\right\} \end{equation}

with $T$ being the amplitude-transmission coefficient, $k_m$ and $k_p$ are the wavenumbers of the wave upstream and downstream of the ramp, respectively. At the points of discontinuity in the bottom slope, namely, $x=0$ and $x=L$, a continuity of both $(a-{1}/{\sigma }) \psi _{0_{x}} + b \psi _{0}$ and $\psi _0$ is required

(3.8)\begin{equation} \left.\begin{gathered} \psi_0 = \psi_0^-,\quad \left(a(h^-, {h^\prime}^+)-\frac{1}{\sigma}\right) \psi_{0_{x}} + b(h^-, {h^\prime}^+,{h^{\prime \prime}}^+) \psi_{0} \\ =\left(a(h^-,0)-\frac{1}{\sigma}\right) \psi_{0_{x}}^- \quad {\rm at}\ x=0,\\ \psi_0 = \psi_0^+,\quad \left(a(h^+, {h^\prime}^-)-\frac{1}{\sigma}\right) \psi_{0_{x}} + b(h^+, {h^\prime}^-,{h^{\prime \prime}}^-) \psi_{0} \\ =\left(a(h^+,0)-\frac{1}{\sigma}\right) \psi_{0_{x}}^+ \quad {\rm at}\ x=L, \end{gathered}\right\} \end{equation}

equating the terms on both sides of the discontinuity points. The $(+)$ and $(-)$ signs indicate taking a limit downstream and upstream of the discontinuity. Normalizing and substituting the streamfunction solution in the flat region into the equations above yields

(3.9)\begin{gather} \psi_{0_x} + \left({\rm i}k_m \frac{a(h^-,0)-\dfrac{1}{\sigma}}{a(h^-, h'^+)-\dfrac{1}{\sigma}} + \frac{b(h^-, h'^+,h''^+)}{a(h^-, h'^+)-\dfrac{1}{\sigma}} \right) \psi_0 = 2{\rm i}\omega A \frac{a(h^-,0)-\dfrac{1}{\sigma}}{a(h^-, h'^+)-\dfrac{1}{\sigma}},\quad x=0, \end{gather}

(3.10)\begin{gather}\psi_{0_x} - \left(ik_p \frac{a(h^+,0)-\dfrac{1}{\sigma}}{a(h^+, h'^-)-\dfrac{1}{\sigma}} + \frac{b(h^+, h'^-,h''^-)}{a(h^+, h'^-)-\dfrac{1}{\sigma}}\right) \psi_0 = 0,\quad x=L. \end{gather}

These are the jump conditions utilized in the numerical solution of the PCMSE.

Figure 8 demonstrates that the reflection coefficient obtained using the PCMSE shows a better agreement with the exact solution than that obtained using the CMSE in the range ${\omega ^2 L}/{g}<0.6$. The MMSE provides the best results in this region. All models are almost indistinguishable in the rest of the region. This result shows yet again the compatibility of the PCMSE to the CMSE for milder slopes and a clearly improved solution for steeper slopes. All in all, the PCMSE solution provides a relatively good agreement in the entire region. It is noted that the left edge of the figure depicts very large bottom slopes which are outside the scope of relevance of these models.

3.4. Bragg resonance – class II

Bragg resonance refers to the interaction of free-surface incident waves with periodic bottom undulations satisfying resonant conditions. This phenomenon plays an important role in the evolution of nearshore waves. The class I Bragg resonance involves interactions of two surface waves and one bottom undulation satisfying the conditions

(3.11)\begin{equation} \left.\begin{gathered} k_1-k_2-k_b = 0,\\ \omega_1 - \omega_2 = 0. \end{gathered}\right\} \end{equation}

The subscript indices refer to the number of the free-surface wave, whereas $b$ refers to the bed. This type of resonance has a significant effect on linear waves. The class II Bragg resonance involves the interaction of two surface waves and two bottom undulations satisfying the conditions

(3.12)\begin{equation} \left.\begin{gathered} k_1-k_2-(k_{b_1} \pm k_{b_2}) = 0,\\ \omega_1 - \omega_2 = 0. \end{gathered}\right\} \end{equation}

This type of resonance has a weaker effect on linear waves than the class I Bragg resonance, but can still have significant effects due to undulating bottom components. The classical MSE, which account for second-order terms, generally provide a good agreement with the class I Bragg resonance. On the other hand, they fail to capture the results of the higher-order conditions of the class II and higher.

The consideration of higher-order bed slopes in the PCMSE is expected to provide a better match to the class II Bragg resonance as compared with the classical MSE. The leading-order term for this class in the bottom slope is the $(\nabla h)^2$ term. The comparison is performed below for small values of the dispersion parameter, $kh$, which is where the PCMSE equation is expected to provide the most accurate results.

Figures 9 and 10 depict the simulation results for the PCMSE, CMSE and MMSE in modelling the class II Bragg resonance reflection coefficient $R$ on a doubly sinusoidal bottom as a function of $\sigma h_0$. The range of $\sigma h_0$ values is $0.02$ to $0.58$ in figure 10 which is within the range where the PCMSE is expected to work well. The range in figure 9 is $0.35$ to $2.15$ which extends well beyond this range and is therefore expected to provide more of a challenge to the models. The simulations of the MSE models are compared with the ‘linear exact’ result, as displayed in figure 3 of Kim et al. (Reference Kim, Ertekin and Bai2010) in figure 9. The first peak located at ${2 k}/{k_{b_1}} \approx 0.3$ in the figure is a result of the class II Bragg resonance. The result of the PCMSE is significantly closer to the linear exact result than the CMSE. The MMSE provides the closest value for this peak. It is noted that Kim et al. (Reference Kim, Ertekin and Bai2010) obtained slightly different results for the CMSE which are closer to the PCMSE results presented here.

Figure 9. Comparison of the values of the reflection coefficient $R$ obtained from simulations of the analytical PCMSE (blue – smaller dashed curve), PCSSE (orange – larger dashed curve), CMSE (red – dot dashed curve) and MMSE (green – dotted curve) along with the result arising from the ‘linear exact’ solution presented by Kim, Ertekin & Bai (Reference Kim, Ertekin and Bai2010) (black – solid curve) taken from figure 3 there. The bed wavenumbers are $k_{b_1} = 2 {\rm \pi}$ cm$^{-1}$ and $k_{b_2} = {4 {\rm \pi}}/{3}$ cm$^{-1}$. The patch length is $L = 12$ cm, the amplitude ratio is ${{\rm \Delta} H}/{H_0} = 0.25$ and $H_0 = 2.5$ cm.

Figure 10. Comparison of the values of the reflection coefficient $R$ obtained from simulations of the analytical PCMSE (blue – smaller dashed curve), PCSSE (orange – larger dashed curve), CMSE (red – dot dashed curve) and MMSE (green – dotted curve) along with the numerical results of Guazzelli, Rey & Belzons (Reference Guazzelli, Rey and Belzons1992) (black – solid curve) taken from figure 2 there. The bed wavenumbers are $k_{b_1} = {{\rm \pi} }/{6}$ cm$^{-1}$ and $k_{b_2} = {{\rm \pi} }/{3}$ cm$^{-1}$. The patch length is $L = 48$ cm, the amplitude ratio ${{\rm \Delta} H}/{H_0} = 0.25$ and $H_0 = 2.5$ cm. (a) An entire domain. (b) Blowup of class II resonance peak.

Figure 10 presents a comparison between the numerical results obtained by Guazzelli et al. (Reference Guazzelli, Rey and Belzons1992) accounting for three evanescent modes on one hand and the results obtained from the MSE on the other hand. The black curve is a reconstruction of the results presented in figure 2 in Guazzelli et al. (Reference Guazzelli, Rey and Belzons1992). The first major peak at $f \approx 1.7$ Hz is a result of the class II Bragg resonance. Figure 10(b) displaying a blowup of this peak shows that the PCMSE-based result is significantly closer to the ‘linear exact’ result than those of the other models.

3.5. Obliquely incident waves

Ehrenmark (Reference Ehrenmark1998) presented an exact solution to obliquely incident waves over a plane beach. The solution is limited, as is the solution of Stoker (Reference Stoker1957), to bottom slopes of ${{\rm \pi} }/{2n}$, where $n$ is an integer. Simulations were performed for a 45 degree bottom slope along with a deep water, 45 degree, oblique incidence angle. Equation (2.41) is utilized for this simulation. The exact solution is compared with the results of numerical simulations for the PCSSE, PCMSE and CMSE in figure 11. A blowup of the wave peak is presented in the bottom part of the figure. Both PC models provide a closer approximation to the peak value than the CMSE when comparing with the exact solution. Closer to the shore, the CMSE provides a slightly better result than the PCSSE with the PCMSE farther off. This difference is diminished among approaching the shore. As an overall observation it seems that the PCSSE provides the most accurate description out of all the simulated MSE models. This simulation provides an initial insight as to the performance of the PC equations for the 3-D case. A more thorough examination is required.

Figure 11. Normalized free-surface values for an oblique incidence simulation for a $45^{\circ }$ bottom plane and a $45^{\circ }$, deep water, incidence angle. (a) The full domain. (b) Blowup of the wave peak. Blue (smaller dashed) curve – PCMSE; orange (larger dashed) curve – PCSSE; red (dot dashed) curve – CMSE; black (solid) curve – exact solution.