2.1. Problem formulation

We consider a classical problem for linearised potential flow (see Stoker Reference Stoker2019) describing the wave field on a beach with constant slope with angle $\alpha$ between the surface of the free water and the beach (bottom of the ocean). In the fluid domain $\varOmega$ (see figure 1), given by

(2.1) \begin{align} \varOmega =\{(r,\theta ): \quad 0\lt r\lt \infty ,\quad \pi \lt \theta \lt \theta _0=\pi +\alpha \leqslant 2\pi \}, \end{align}

the velocity potential $\phi (r,z)$ satisfies the Laplace equation

(2.2) \begin{equation} \Delta \phi = 0. \end{equation}

Due to the problem geometry, it is convenient to consider it in polar coordinates

(2.3) \begin{align} x=r\cos \theta ,\quad z=r\sin \theta , \end{align}

with the origin at the corner where the fluid surface and the beach meet. Along the beach, the Neumann boundary condition of no flow applies, so that

(2.4) \begin{equation} \frac {\partial }{\partial n}\phi (r,\theta )\big |_{\theta =\theta _0} = 0, \end{equation}

where $n$ is the outer normal to the boundary. In accordance with the linearised theory, at the water surface, we define the surface displacement by $\zeta (t,x)$ . The kinematic condition

(2.5) \begin{align} \partial _z\phi = \partial _t \zeta , \quad z=0, \end{align}

and the dynamic condition (Mei Reference Mei1989)

(2.6) \begin{align} \partial _t\phi + g \zeta =0,\quad z=0, \end{align}

where $g$ is the gravitational acceleration. In the frequency domain, assuming all functions are proportional to $\mathrm{e}^{-\mathrm{i} \omega t}$ , these conditions become

(2.7) \begin{equation} \partial _z \phi =k\phi ,\quad z=0, \end{equation}

where $k = \omega ^2/g$ and $\omega$ is the angular frequency.

This problem has been addressed using the complex analysis in Lewy (Reference Lewy1946), Stoker (Reference Stoker1947), Isaacson (Reference Isaacson1950) and Peters (Reference Peters1950), while Ehrenmark (Reference Ehrenmark1987, Reference Ehrenmark1989) found a closed form solution representation using Mellin transform. In the first two papers, exact solutions have been constructed for domains with specific angles, respectively, $\theta _0=\pi +\pi / (2m)$ and $\theta _0=\pi +\pi p/(2m)$ ( $p\lt 2m$ , where $p$ and $m$ are natural numbers), while the latter extended the method to consider an arbitrary domain $\pi \lt \theta _0\leqslant 2\pi$ . Note that, in their analysis, they measure the angle differently: $0\lt \theta \lt \alpha \leqslant \pi$ ( $\theta _0=\alpha +\pi$ ).

Figure 1.

Schematic diagram showing the equations of motion and boundary conditions.

It was shown that this linear homogeneous problem has a non-trivial solution that describes sinusoidal waves along the water surface at infinity ( $r\to \infty$ ) and it is bounded as $r\to 0$ (see figure 1) (we note that the classic papers mentioned above, in addition to the regular solution considered in this paper, also discussed singular solutions; we will consider those in a separate paper)

(2.8) \begin{equation} \phi (r,\theta )=\phi _0 +\phi _1 r^{\gamma _0}+O\left (r^{\gamma _0+\delta }\right )\!,\quad r\to 0, \quad \pi \lt \theta \leqslant \theta _0\leqslant 2\pi , \end{equation}

Here, $\phi _{\!j}=\phi _{\!j}(\theta _0)$ and $0\lt \gamma _0(\theta _0)\leqslant 1$ are to be determined, while $\delta \gt 0$ is a small number.

To ensure the uniqueness of the solution, we require the following asymptotic estimate is valid:

(2.9) \begin{equation} \phi (x,z) = e^{i k x}e^{k z} + R\, e^{-i k x}e^{k z}+o(1),\quad x\to -\infty . \end{equation}

Here, $R$ is the complex reflection coefficient that in this case represents a shift in the mode between the incident and the reflected waves (which are the first and second terms in (2.9), respectively). For the regular solution, energy conservation gives $|R|=1$ . Moreover, for those analytical solutions presented in Lewy (Reference Lewy1946), the remainder (the $o(1)$ term in (2.9)) exponentially decays at infinity and represents evanescent waves. However, when considering the obtuse beach angle ( $\pi /2\lt \alpha \lt \pi$ ) or dock problem ( $\alpha =\pi$ ) it is considered that the decay is slower (algebraic).

We do not require any of these facts (knowledge of the reflection coefficient, $R$ , or the rate of decay) and will receive the proper results in due course of our analysis.

Remark 1. Note that, for the vertical slope ( $\theta _0=3\pi /2$ , or equivalently, $\alpha =\pi /2$ ) the $\text{regular}$ solution is trivial and represents a well-known standing wave

(2.10) \begin{equation} \phi (x,z) = e^{i k x}e^{k z} + e^{-i k x}e^{k z}=2e^{k z}\cos kx, \end{equation}

thus, $R=1$ in this case.

In the next subsection, we reformulate the above problem using the Mellin transform, paying particular attention to the properties of analytical function that represents the Mellin image of the potential $\phi (r,\pi )$ along the surface of the free water. Those properties are natural and are crucial in understanding our further transformations.

2.2. Transformation of the problem using Mellin transform

Following Ehrenmark (Reference Ehrenmark1987), we apply the Mellin transform, defined as

(2.11) \begin{align} \mathcal{M}\{\phi \}(s,\theta ) \equiv \tilde \phi (s,\theta )= \int _0^\infty r^{s-1} \phi (r, \theta ) \, \text{d}r \end{align}

to the governing (2.2) and the boundary conditions (2.4) and (2.7). Here $s$ is the Mellin transform variable, a complex number belonging to a domain to be discussed later. This leads to the trivial ordinary differential equation

(2.12) \begin{equation} \tilde \phi _{\theta \theta } +s^2 \tilde \phi =0. \end{equation}

According to estimates (2.8) and (2.9), the solution sought does not belong to a classic functional space of this transform. Thus, we should operate in a space of distribution (Zemanian Reference Zemanian1987). This is not difficult in this case due to the smoothness of the solution inside the domain but requires the respective justifications.

To avoid this route, we consider a regularised solution to the problem with a slightly deformed asymptotic estimate at infinity

(2.13) \begin{equation} \phi (r,\pi )=e^{-ik r}e^{-\epsilon r} + R\, e^{ik r}e^{-\epsilon r}+O\big (r^{-\gamma _\infty }\big ),\quad r\to \infty , \end{equation}

where $\epsilon \gt 0$ can later be taken to zero to allow both the incident and a reflected wave to be present. Again, the parameters $R$ and $\gamma _\infty \gt 0$ are considered unknown.

As a result, we expect that the Mellin image, $\tilde \phi (s,\theta )$ , is an analytic function in the strip $0\lt \mathrm{Re} s\lt \gamma _\infty$ for any $\pi \lt \theta \leqslant \theta _0$ and $\pi \lt \theta _0\leqslant 2\pi$ .

Moreover, a priori estimates (2.8) and (2.13) suggest that the function can possess simple poles

(2.14) \begin{equation} \tilde \phi (s,\theta )\sim \frac {\phi _0}{s}+O(1), \quad s\to 0, \end{equation}

while

(2.15) \begin{equation} \tilde \phi (s,\theta )\sim \frac {\phi _1}{s+\gamma _0}+O(1), \quad s\to -\gamma _0, \end{equation}

and the function $\tilde \phi (s,\theta )$ can be analytically extended into the strip

(2.16) \begin{equation} -\gamma _0-\delta \lt \mathrm{Re} s\lt \gamma _\infty , \end{equation}

where $\delta \gt 0$ .

Remark 2. We expect the following for the unknown parameter that defines the estimate of the solution at zero and infinity:

(2.17) \begin{equation} \begin{matrix} \gamma _0=1,\quad \gamma _\infty =\infty , &\mathrm{when} \quad 0\lt \alpha \lt \pi /2, \\[3mm] 0\lt \gamma _0(\alpha )\lt 1, \quad \gamma _\infty (\alpha )\gt 1, &\mathrm{when} \quad \pi /2\lt \alpha \lt \pi . \end{matrix} \end{equation}

In what follows, we refer to the unknown strip ( 2.16 ), specifying the parameters when appropriate.

Finally, we observe that the function $\tilde \phi (it,\theta )$ should decay faster than any inverse polynomial as $t\to \pm \infty$ , due to the smoothness of the original solution.

Using the boundary condition (2.4) we can conclude that

(2.18) \begin{equation} \tilde \phi (s,\theta )=A(s)\cos (\theta -\theta _0) s, \end{equation}

or

(2.19) \begin{equation} \phi (r,\theta )=\frac {1}{2\pi i}\int _{-i\infty +\eta }^{i\infty +\eta } A(s) \cos \big (s(\theta -\theta _0)\big )r^{-s}ds,\quad 0\lt \eta \lt \gamma _\infty . \end{equation}

An equivalent representation of (2.19) is

(2.20) \begin{equation} \phi (r,\theta )=\frac {1}{2\pi i}\int _{-i\infty +\eta }^{i\infty +\eta } \tilde \phi (s,\pi ) \frac {\cos s(\theta -\theta _0)}{\cos s(\pi -\theta _0)}r^{-s}ds, \end{equation}

demonstrating that the integration path can be shifted to the left, revealing the asymptotic estimate (2.8).

In the next section, we formulate a functional equation (2.33) for an auxiliary function $w(r)$ (its Mellin image is an analytical function $\tilde w(s)$ ) that preserves the properties of the potential along the free water surface ( $\theta =\pi$ ) and decays at infinity (thus the incident and reflected waves are participated in the analysis separately in their analytical forms).

2.3. Transforming the original problem to a functional equation

Let us introduce a representation of the solution on the surface $z=0$ in the form

(2.21) \begin{equation} \phi (r,\pi )=U_\epsilon (r)+w(r),\quad U_\epsilon (r)=e^{-ik r}e^{-\epsilon r} + R\, e^{ik r}e^{-\epsilon r}, \end{equation}

where $U_\epsilon$ represents the incident and reflected wave contribution to the surface potential. We can write (compare with (2.8) and (2.9)) we can write

(2.22) \begin{align} w(r)=w_0+ w_1r^{\gamma _1}+O(r^{\gamma _1+\delta }),\quad r\to 0, \\[-12pt]\nonumber \end{align}

(2.23) \begin{align} w(r)= w_\infty r^{-\gamma _\infty }+O(r^{-\gamma _\infty -\delta }),\quad r\to \infty , \end{align}

where

(2.24) \begin{equation} \phi _0=1+R+w_0,\quad \gamma _1=\min \{1,\gamma _0\}\gt 0, \end{equation}

and both constants $R$ and $w_0$ are unknown at this stage and will be determined later.

Note that determining the function $w(r)$ is enough to solve the original problem due to (2.21) and the representation (2.20).

We emphasise that both functions $\phi$ and $w$ actually depend on the parameter $\epsilon$ (that is, $\phi = \phi _\epsilon$ and $w = w_\epsilon$ ), while the solution to the original problem will be obtained in the limit $\epsilon \to 0$ . For the sake of clarity, we will omit the additional subscripts and simply indicate where the limit is taken to recover the solution.

Next, we compute

(2.25) \begin{equation} \tilde \phi (s,\pi )=\tilde w(s)+\tilde U_\epsilon (s). \end{equation}

Moreover (compare with (2.22)),

(2.26) \begin{equation} \tilde w(s)=\frac {w_0}{s}+O(1),\quad s\to 0, \end{equation}

(2.27) \begin{equation} \tilde w(s)=\frac {w_1}{s+\gamma _0}+O(1),\quad s\to -\gamma _0, \end{equation}

and there are no other poles of the function in the strip $-\gamma _1-\delta \lt \mathrm{Re} s\lt \gamma _\infty$ .

Remark 3. In the following, we assume

(2.28) \begin{equation} \gamma _1+\gamma _\infty +\delta \gt 1, \end{equation}

that will be justified in due course.

As a result, functions $s\tilde w(s)$ and $\tilde w(s+1)$ are analytical in the common strip

(2.29) \begin{equation} -\gamma _1-\delta \lt \mathrm{Re} s\lt \gamma _\infty -1, \end{equation}

having no poles except at those points mentioned in ( 2.26 ) and ( 2.27 ), and they both decay at infinity as

(2.30) \begin{equation} \tilde w(s)=O(|s|^{-\beta }),\quad |\mathrm{Im} s|\to \infty . \end{equation}

Here, $\beta \gt 0$ is an arbitrary constant.

The other term in (2.25) can be computed exactly to read

(2.31) \begin{equation} \tilde U_\epsilon (s)= \frac {\Gamma (s)}{(\epsilon +ik)^{s}}+ \frac {\Gamma (s)}{(\epsilon -ik)^{s}}R. \end{equation}

Note that the latter is analytical in the strip $0\lt \mathrm{Re} s \lt b$ ( $b\gt 0$ is an arbitrary value) and $\tilde U_\epsilon (s)$ has a simple pole at zero

(2.32) \begin{equation} \tilde U_\epsilon (s)\sim \frac {1+R}{s},\quad s\to 0, \end{equation}

while $\tilde U_\epsilon (s)$ has a simple pole at the point $s=-1$ (related to the behaviour of the gamma function). Both functions $\tilde U_\epsilon (s)$ and $\tilde U_\epsilon (s+1)$ decay at infinity faster than any inverse polynomial of $s$ ( $|\mathrm{Im} s|\to \infty$ ) provided $\epsilon \gt 0$ .

Applying the Mellin transform to the remaining unused condition on the free water surface (2.7), we obtain a functional equation

(2.33) \begin{equation} s\tilde w(s)\tan \big (s(\pi -\theta _0)\big )-k\tilde w(s+1)=\tilde f_\epsilon (s), \end{equation}

valid in the strip (2.29) where

(2.34) \begin{equation} \tilde f_\epsilon (s)=- s\tilde U_\epsilon (s)\tan (s(\pi -\theta _0))+k\tilde U_\epsilon (s+1). \end{equation}

Let us now pass to the limit $\epsilon \to 0$

(2.35) \begin{equation} \tilde U_0(s)= \lim _{\epsilon \to 0}\tilde U_\epsilon (s)= \frac {\Gamma (s)}{k^{s}} \left ((1+R)\cos \frac {\pi s}{2}-i(1-R)\sin \frac {\pi s}{2}\right )\!, \end{equation}

while

(2.36) \begin{equation} \tilde U_0(s+1)= \frac {\Gamma (1+s)}{k^{s+1}} \left (-(1+R)\sin \frac {\pi s}{2}-i(1-R)\cos \frac {\pi s}{2}\right )\!. \end{equation}

This allows us to compute a limiting function from (2.34)

(2.37) \begin{equation} \tilde f_0(s)=- \dfrac {\Gamma (s+1)}{k^{s}} \left ( (1+R) \dfrac {\sin s\left(\dfrac {3\pi }{2}-\theta _0\right)}{\cos s(\pi -\theta _0)}+i (1-R)\dfrac {\cos s\left(\dfrac {3\pi }{2}-\theta _0\right)}{\cos s(\pi -\theta _0)} \right )\!. \end{equation}

Note that $\tilde f_0(s)$ is an analytical meromorphic function in the strip containing the origin $s=0$

(2.38) \begin{align} -\min \left \{1,\frac {\pi }{2(\theta _0-\pi )}\right \}\lt \mathrm{Re} s \lt \frac {\pi }{2(\theta _0-\pi )}, \end{align}

and can be considered in any wider strip where it has poles related to the gamma function ( $s=-1,-2,\ldots$ ) and the trigonometric functions ( $s=(({\pi +2\pi j})/({2(\theta _0-\pi )}))$ , $j=\pm 1,\pm 2,\ldots$ ). Moreover, for any $\theta _0\gt \pi$ , $\tilde f_0(s)$ exponentially decays at infinity within the strip. In this way, we have justified our notation. Indeed, there exists a function $f_0(r)$ such that $\tilde f_0(s)={\mathcal M}[f_0](s)$ .

Rewriting (2.33) in the form

(2.39) \begin{equation} s \tilde w(s)=\frac {k\tilde w(s+1)+\tilde f_0(s)}{\tan \big (s(\pi -\theta _0)\big )}, \end{equation}

analysing it as well as (2.33) we conclude in the worst-case scenario that

(2.40) \begin{equation} \gamma_0=\min\left\{1,\frac{\pi}{2(\theta_0-\pi)}\right\}\!,\quad \gamma _\infty =1+\frac {\pi }{2(\theta _0-\pi )}, \end{equation}

thus the assumption (2.28) is always satisfied. Note that the strip (2.29) can be wider (for example, if some of the constants $w_1$ or $w_\infty$ vanishes). Moreover, we can always analytically extend the equation in any finite strip where the equation can have isolated poles. Finally, by analysing (2.39) near the point $s=0$ , we conclude

(2.41) \begin{equation} k\tilde w(s+1)+ k\tilde U_0(s+1)\sim s(\pi -\theta _0)w_0 ,\quad s\to 0. \end{equation}

Note that

(2.42) \begin{equation} \tilde w(1)= -\tilde U_0(1)=\frac {i(1-R)}{k}, \end{equation}

represents the solvability condition for the Laplace equation in the wedge domain with the Neumann boundary condition on the sides on the edges (see figure 1).

From this point on, the functions $w(r)$ and $\tilde w(s)$ are related to the original solution obtained in the limit $\epsilon \to 0$ .

While the Mellin-transform framework employed up to this point can naturally share elements with the analyses presented in Ehrenmark’s papers, the subsequent development introduces a distinct formulation that diverges from his methodology. That is, we are now going to further transform (2.33) to produce another functional equation (2.48). The second equation differs from the first one in that all the coefficients in front of the unknown functions have a finite, non-zero limit as their argument tends to infinity along the imaginary axis in the complex plane. This, in turn, allows us to transform the latter to an integral equation with a kernel that admits a good behaviour along the real axis.

2.4. Transformation of functional equation (2.33)

To achieve the goal, we use a similar line of reasoning as in Mishuris (Reference Mishuris2001) by looking for a solution of (2.33) in terms of a new function $p(r)$

(2.43) \begin{equation} \tilde w(s+1)=-\dfrac {1}{k^{s+1}}\Gamma (s+1)\sin \dfrac {\pi s}{2} \left (\tilde p(s+1)-\dfrac {A}{\sin \dfrac {\pi s}{2} }\right )\!, \end{equation}

(2.44) \begin{equation} \tilde w(s)=\dfrac {1}{k^s}\Gamma (s)\cos \dfrac {\pi s}{2}\left (\tilde p(s)+\dfrac {A}{\cos \dfrac {\pi s}{2} }\right )\!, \end{equation}

where function $\tilde p(s)$ is the Mellin transform of the unknown function $p(r)$ which is analytic on the same strip with some change in behaviour near the points $s=0,\pm 1$ , while the constant $A$ will be defined later.

Remark 5. Note that, for any real $t$ ,

(2.45) \begin{align} \left |\Gamma (it)\right |=\sqrt{\frac{\pi}{t\sinh \pi t}},\quad \end{align}

and thus the function $\tilde p(s)$ introduced in ( 2.44 ) behaves at infinity in a manner similar to the original function $\tilde w(s)$

(2.46) \begin{equation} \frac{\tilde{w}(it)}{\tilde{p}(it)}=O(t^{-1/2}),\quad t\to\pm\infty \end{equation}

thus both of them should simultaneously decay at infinity or

(2.47) \begin{equation} \lim _{t\to \pm \infty } {\tilde p(it)}=0 \end{equation}

(compare this with the assumption ( 2.30 )). Lastly, the function $\tilde p(s)$ , in contrast to $\tilde w(s)$ , does not have a pole at the point $s=0$ .

Substituting (2.43)–(2.44) into (2.33) we arrive at a new functional equation

(2.48) \begin{equation} \tan \big (s(\theta _0-\pi )\big )\tilde p(s) - \tan \frac {\pi s}{2} \tilde p(s+1)=\tilde f_1(s), \end{equation}

in the strip

(2.49) \begin{equation} -1\lt \ \mathrm{Re} s \lt \frac {\pi }{2(\theta _0-\pi )}, \end{equation}

where the functions $\tilde p(s)$ and $\tilde p(s+1)$ do not have any poles. The function $\tilde f_1(s)$ is given by

(2.50) \begin{equation} \tilde f_1(s)=-\frac {1}{\cos \dfrac {\pi s}{2}}\left ( \dfrac {k^s}{\Gamma (s+1)}\tilde f_0(s)+A\left (\tan \big (s(\theta _0-\pi )\big )+1\right )\right )\!, \end{equation}

the constant $A$ can now be evaluated so as to ensure that $\tilde p(s)$ is analytic at $s=1$ (compare (2.26), (2.41), (2.36) and (2.44))

(2.51) \begin{equation} A=i(1-R),\quad \tilde p(0)+A=w_0. \end{equation}

To summarise the rationale behind the substitution (2.43), we note that the new unknown function $\tilde p(s)$ (i) exhibits the same asymptotic behaviour at infinity within the strip of analyticity as $\tilde w(s)$ , (ii) does not possess poles in a wider strip (does not have the pole at $s=0$ ) and (iii) appears in the functional equation (2.33) without algebraically growing coefficients in front of the unknowns $\tilde p(s)$ and $\tilde p(s+1)$ . This allows us to apply the inverse Mellin transform to the equation using a Fubini-type theorem and avoids the need to work in distributional spaces. Finally, we have eliminated the gamma function from the right-hand side of the functional equation (2.33) (compare (2.37) and (2.50)).

We observe that the right hand side of (2.48) can be split into the sum

(2.52) \begin{equation} \tilde f_1(s)=(1+R)\tilde f_1^{(1)}(s)+i(1-R)\tilde f_1^{(2)}(s), \end{equation}

where

(2.53) \begin{align} \tilde f_1^{(1)}(s) &= \dfrac {\sin s \left(\dfrac {3\pi }{2}-\theta _0\right)}{\cos s(\pi -\theta _0)\cos \dfrac {\pi s}{2}},\nonumber \\ \tilde f_1^{(2)}(s) &= \dfrac {\cos s \left(\dfrac {3\pi }{2}-\theta _0\right)}{\cos s(\pi -\theta _0)\cos \dfrac {\pi s}{2}}- \dfrac {\tan (s(\theta _0-\pi ))+1}{\cos \dfrac {\pi s}{2}}. \end{align}

This representation allows us to search for the solution of (2.48) in a form similar to (2.52)

(2.54) \begin{equation} \tilde p(s)=(1+R)\tilde p^{(1)}(s)+i(1-R)\tilde p^{(2)}(s), \end{equation}

where each of the functions with the subscript $j=1,2$ satisfy the respective equation

(2.55) \begin{equation} \tan s (\theta _0-\pi )\tilde p^{(j)}(s)- \tan \frac {\pi s}{2} \tilde p^{(j)}(s+1)=\tilde f^{(j)}_1(s). \end{equation}

Note that functions $f_1^{(j)}(s)$ can be estimated at infinity inside the analyticity strip

(2.56) \begin{equation} \tilde f_1^{(j)}(s)=O(e^{q|s|}), \quad \mathrm{Im} s\to \pm \infty , \end{equation}

with $q=|3\pi /2-\theta _0|-\theta _0+\pi /2\leqslant -\pi$ .

Moreover,

(2.57) \begin{align} \tilde f^{(2)}_1(s)=-s(\theta _0-\pi )+O(s^2),\quad \tilde f^{(1)}_1(s)=s \left (\frac {3\pi }{2}-\theta _0\right )+O(s^3), \quad s\to 0. \end{align}

Thus, since $\tilde p^{(j)}(s)$ is analytic near the point $s=0$ , and the same happens at $s=1$ , the left-hand side of (2.55) vanishes as $s\to 0$ .

As a result of this subsection, we have arrived at two independent functional equations (2.55) that are essentially the same (having only different right-hand sides). When we further transform them into integral equations with fixed-point singularities at the ends of the integration path, we will observe that the solution may possess an additional unwanted property. Namely, they might have a singularity $(r-1)^{-1}$ , as $r\to 1$ and thus their analysis in the classic functional spaces natural for such equations cannot be performed. We try to highlight the problem differently. Note that the auxiliary function $\tilde p(s)$ has to satisfy the additional condition (2.47), while the functions $\tilde p^{(j)}(s)$ do not necessarily need to possess such properties. For example, it is easy to observe that any of the assumptions

(2.58) \begin{equation} \tilde p^{(j)}(it)\to C_{\!j},\quad \mbox{or}\quad \tilde p^{(j)}(it)\to C_{\!j} \tan (it),\quad t\to \pm \infty \end{equation}

do not contradict (2.55) and this fact should be taken into account.

2.5. Further transformation of functional equations (2.55)

To address this fact, we represent the unknown solutions as follows:

(2.59) \begin{equation} p^{(j)}(r)= \frac {2d_{\!j}}{\pi }\left (\frac {r}{r^2-1}-\frac {r}{r^2+1}\right )+p^{(j)}_*(r), \end{equation}

or

(2.60) \begin{align} \tilde p^{(j)}(s)&= d_{\!j}\left (\tan \dfrac {\pi s}{2}-\dfrac {1}{\cos \dfrac {\pi s}{2}}\right )+\tilde p^{(j)}_*(s),\\[3pt] \tilde p^{(j)}(s+1)&= -d_{\!j}\left (\cot \dfrac {\pi s}{2}-\dfrac {1}{\sin \dfrac {\pi s}{2}}\right )+\tilde p^{(j)}_*(s+1), \nonumber \end{align}

where functions $p^{(j)}_*(r)$ are regular near the point $r=1$ or, equivalently,

(2.61) \begin{equation} \tilde p^{(j)}_*(s)\to 0, \quad |\mathrm{Im} s|\to \infty . \end{equation}

Thus, instead of two functions $\tilde p^{(j)}(s)$ , we need to find two constants $d_{\!j}$ and functions $\tilde p^{(j)}_*(s)$ . Note that $\tilde p^{(j)}(s)$ from (2.60) generally do not decay at infinity. To guarantee the decay of the function $\tilde p(s)$ (compare (2.47), (2.54) and (2.61)) we put

(2.62) \begin{equation} (1+R)d_1+i(1-R)d_2=0, \end{equation}

that, in turn, defines the reflection coefficient $R$ . We would like to underline here that nothing at this stage shows that condition $|R|=1$ is valid, since the values of the (complex) constants $d_{\!j}$ are unknown.

The final solution $\tilde p(s)$ takes the form (note that $p^{(j)}_*(s)$ and $p^{(j)}(s)$ from (2.54) are different)

(2.63) \begin{equation} \tilde p(s)=(1+R)\tilde p^{(1)}_*(s)+i(1-R)\tilde p^{(2)}_*(s). \end{equation}

As a result, the functional equations (2.55), written in terms of the functions $\tilde p^{(j)}_*(s)$ , take similar forms to the other right-hand sides and while their solutions are bounded as $s\to \infty$ inside the analyticity strip (2.49)

(2.64) \begin{equation} \tan s (\theta _0-\pi )\tilde p^{(j)}_*(s)- \tan \frac {\pi s}{2} \tilde p_*^{(j)}(s+1)=\tilde f^{(j)}_\clubsuit (s),\quad j=1,2, \end{equation}

where

(2.65) \begin{align} \tilde f^{(j)}_\clubsuit (s)=\tilde f^{(j)}_1(s)-d_{\!j} \tilde f^{(2)}_1(s). \end{align}

After this preparatory work, we are ready to evaluate the integral equations that possess enough flexibility (additional parameters) to eliminate that possible singular behaviour at the point $r=1$ .

2.6. Transformation of functional equations (2.64) to integral equations

We now divide (2.64) by $\tan ({\pi s}/{2})$ and observe that

(2.66) \begin{align} \cot \frac {\pi s}{2} \tan s (\theta _0-\pi )=1-\tilde n_0(s), \end{align}

where

(2.67) \begin{equation} \tilde n_0(s)=\frac {\sin s\left(\dfrac {3\pi }{2}-\theta _0\right)}{\sin \dfrac {\pi s}{2} \cos s(\theta _0-\pi )}. \end{equation}

This allows us, by applying the inverse Mellin transform, to rewrite equations (2.64) in the form

(2.68) \begin{equation} (1-r) p^{(j)}_*(r)-{\mathcal N}\,[p^{(j)}_*](r)=f_{*}^{(j)}(r)\equiv \frac {1}{2\pi i} \int _{-i\infty +\eta }^{i\infty +\eta }\frac {\tilde f_{\clubsuit }^{(j)}(s)r^{-s}\,\text{d}s}{\tan (\pi s/2)}, \end{equation}

where $\eta \gt 0$ is an arbitrary value such that $\eta \lt \min \{\gamma _1,\gamma _\infty -1\}$ , where $\gamma _1$ and $\gamma _\infty$ were evaluated in (2.40). Here, the integral operator $\mathcal N$ is

(2.69) \begin{equation} {\mathcal N}\,[u](r)= \int _{0}^{\infty }n_0(t/r)u(t)\frac {\text{d}t}{t}, \end{equation}

and its kernel, $n_0(r)$ , defined in (2.67) can be simplified for computations

(2.70) \begin{equation} n_0(r)=\frac {1}{\pi } \int _0^\infty \frac {\sinh (3\pi /2-\theta _0)\xi \cos (\xi \log r)}{\sinh (\pi \xi /2) \cosh (\theta _0-\pi )\xi }\text{d}\xi . \end{equation}

To our knowledge, there are no closed analytical representations for this function for an arbitrary value of $\theta$ (see Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2014).

Since the unknown function $p^{(j)}_*(r)$ should be smooth along the positive part of the real axis, we additionally require the following two conditions to be valid to eliminate the possible singularity at the point $r=1$ :

(2.71) \begin{equation} {\mathcal N}\,[p^{(j)}_*](1)+f_{*}^{(j)}(1)=0. \end{equation}

These, in turn, define the two still unknown constants $d_{\!j}$ .

The right-hand sides of (2.68) can be further evaluated

(2.72) \begin{align} f_{*}^{(1)}(r)\equiv \frac {1}{2\pi i} \int _{-i\infty +\eta }^{i\infty +\eta } \left (\tilde f^{(1)}_1(s)- d_1\tilde f^{(2)}_1\right )\frac {r^{-s}\text{d}s}{\tan (\pi s/2)}=n_0(r)-d_1 n_d(r), \end{align}

where

(2.73) \begin{equation} n_d(r)= \dfrac {1}{2\pi i} \int _{-i\infty }^{i\infty } \left ( \dfrac {1}{\sin \dfrac {\pi s}{2}} \left ( \dfrac {\cos s\left(\dfrac {3\pi }{2}-\theta _0\right)}{\cos s(\theta _0-\pi )}-1\right )- \dfrac {\tan s(\theta _0-\pi )}{\sin \dfrac {\pi s}{2}}\right )r^{-s}ds, \end{equation}

while

(2.74) \begin{align} f_{*}^{(2)}(r)\equiv \frac {1}{2\pi i} \int _{-i\infty +\eta }^{i\infty +\eta }\big (\tilde f^{(2)}_1(s)- d_2f^{(2)}_1(s)\big )\frac {r^{-s}\text{d}s}{\tan (\pi s/2)} =(1-d_2)n_d(r). \end{align}

Note that for $\alpha \ne \pi /2$ , any natural $j=1,2,3,\ldots$ , and $l=0,d$

(2.75) \begin{equation} t^j\frac {\text{d}^j}{\text{d}t^j}n_l(t)=O(t^{\beta _*}),\quad t\to 0,\quad t^j\frac {\text{d}^j}{\text{d}t^j}n_l(t)=O(t^{-\beta _*}),\quad t\to \infty , \end{equation}

where

(2.76) \begin{equation} \beta _*=\min \left \{1,\frac {\pi }{2\alpha }\right \}=\min \left \{1,\frac {\pi }{2(\theta _0-\pi )}\right \}\!. \end{equation}

Then, (2.68) read (introducing the operator $\mathcal B$ )

(2.77) \begin{equation} {\mathcal B}[p^{(j)}_*](r)\equiv (1-r) p^{(j)}_*(r)-{\mathcal N}\,[p^{(j)}_*](r)=n_{\!j}(r)-d_{\!j}n_d(r),\quad j=1,2, \end{equation}

where $n_1(r)=n_0(r)$ , $n_2(r)=n_d(r)$ , while the respective conditions (2.71) read

(2.78) \begin{equation} {\mathcal N}\,[p^{(j)}_*](1)+n_{\!j}(1)=d_{\!j}n_d(1). \end{equation}

Note that since

(2.79) \begin{align} n_d(1)=-\frac {1}{\pi } \int _{0}^{\infty } \frac {\tanh t(\theta _0-\pi )}{\sinh \dfrac {\pi t}{2}}\text{d}t\lt 0, \end{align}

conditions (2.78) make sense. However, one can expect that for small values of the angle $\alpha =\theta _0-\pi$ , computations would be challenging (not a surprise, as here the equation of shallow water should take over).

Statement 1. For almost all $0\lt\alpha\le\pi$ , each of equation (2.77), with the corresponding condition (2.78), has a unique smooth solution on the open positive half-axis that satisfies the following estimate at the ends of the interval:

(2.80) \begin{equation} p^{(j)}_*(r)\sim p_0^{(j)} r^{\beta _0},\quad r\to 0,\quad p^{(j)}_*(r)\sim p_\infty ^{(j)} r^{-\beta _\infty },\quad r\to \infty . \end{equation}

The exceptional angles are discussed in Remarks 9 and 10.

The idea of the proof. Details on the analysis of integral equations with fixed-point singularities of that kind can be found in Mishuris (Reference Mishuris1995b ). Moreover, the equation can be solved numerically by the respective projection methods (Gokhberg & Feldman Reference Gokhberg and Feldman1974).

Although the behaviour at infinity is easy to justify since $\beta _\infty =1+\beta _*$ (see (2.76)), the other estimate is more difficult as it also involves analysis of the symbol of operator $\mathcal B$ (see Mishuris Reference Mishuris1995b ). Indeed, the behaviour at the zero point is defined by the zeros of the function $1-\tilde n_0(s)$ that is the symbol of the integral operator $\mathcal B$ defined on the left side of (2.77). Note that, since $1-n_0(0)=2(\theta _0-\pi )/\pi \gt 0$ (compare with (2.67)), it follows that $1-\tilde n_0(it)\gt 0$ for any real $t$ . The zero of the symbol possessing the minimal real part inside the analyticity strip is real, and the exponent in the asymptotic estimate at the zero point, $\beta _0\gt 0$ , is defined from the equation: $ \sin \beta _0(\theta _0-\pi )\cos (\pi \beta _0/2)=0,$ or $\beta _0=\min \{1,\pi /(\theta _0-\pi )\}$ . On the other hand, the right-hand side of (2.77) (and the operator involved) can be estimated, as (2.75) defines the leading-order term that explains (2.80). As a result, we conclude that $\beta _0=\beta _*$ . In other words, the functions $\tilde p_*^{(j)}(s)$ are analytical at are analytical in the strip

(2.81) \begin{equation} -\min\{\beta_0,\beta_*\}\lt\Re s\lt \min\{1+\beta_*,\beta_\infty\}. \end{equation}

Remark 7. Note that both the parameter $d_1$ and the respective smooth solution $p^{(1)}(r,d_1)$ of the remaining integral equation are real. The reflection coefficient should then be found from ( 2.62 )

(2.82) \begin{equation} R=-\frac {d_1+id_2}{d_1-id_2}=\frac {1-id_1}{1+id_1}, \end{equation}

noting that $d_2=1$ . This shows that the reflection coefficient indeed satisfies the condition $|R|=1$ . Moreover, the phase shift $\nu$ is related to the reflection coefficient

(2.83) \begin{equation} R=e^{2i\nu },\quad \tan \nu =-d_1. \end{equation}

Now, when we have proved the solvability of the integral equations, we can reconstruct the solution of the original problem described in § 2.1.

2.7. Reconstruction of the original solution

Taking into account (2.35), (2.44) and (2.51) we have

(2.84) \begin{align} \tilde \phi (s,\pi )&=\tilde U_0(s)+ \frac {1}{k^s}\Gamma (s)\cos \frac {\pi s}{2}\left (\tilde p(s)+\frac {A}{\cos \dfrac {\pi s}{2} }\right ) \nonumber \\ &= \frac {\Gamma (s)}{k^{s}} \left ((1+R)\cos \frac {\pi s}{2}-(1-R)i\sin \frac {\pi s}{2}+ \cos \frac {\pi s}{2}\tilde p(s)+i(1-R)\right )\!. \end{align}

We can now observe that for any $k\gt 0$

(2.85) \begin{align} \int _0^\infty x^{s-1} \sin kx \text{d}x & =k^{-s}\Gamma (s)\sin \frac {\pi s}{2},\quad -1 \lt \mathrm{Re} s \lt 1, \nonumber \\ \int _0^\infty x^{s-1} \cos kx \text{d}x & =k^{-s}\Gamma (s)\cos \frac {\pi s}{2},\quad 0 \lt \mathrm{Re} s \lt 1. \end{align}

Applying the inverse Melling transform to (2.20) (for $\theta =\pi$ ) and any $0\lt \eta \lt \beta _*$ to (2.84), while using the convolution theorem for the Melin transform (see, for example, Zemanian (Reference Zemanian1987)), allows us to reconstruct the potential on the water surface. As discussed in (2.63), $p(r)$ is represented as a sum of two terms; therefore,

(2.86) \begin{equation} \phi (r,\pi )=(1+R)\big (\cos kr+T(kr)\big )-i(1-R)\big (\sin kr- e^{-kr}\big ). \end{equation}

Here we define ( $\eta \gt 0$ )

(2.87) \begin{equation} T(y)=\frac {1}{2\pi i} \int _{-i\infty +\eta }^{i\infty +\eta } \Gamma (s)\cos \frac {\pi s}{2}\tilde p^{(1)}_*(s)y^{-s}\text{d}s= \int _0^\infty \cos \left(\frac {y}{t}\right) p^{(1)}_*(t)\,\frac {\text{d}t}{t}. \end{equation}

Note that solution (2.86) can be written in the form

(2.88) \begin{equation} \phi (r,\pi )=2e^{i\nu }\big ( \cos (kr+\nu )+T(kr) \cos \nu +e^{-kr}\sin \nu \big ). \end{equation}

The value of the potential at the vertex is computed as follows:

(2.89) \begin{equation} \phi (0,\pi )= 2e^{i\nu }\big ( \big (1+\tilde p_*^{(1)}(0)\big )\cos \nu +\sin \nu \big ), \end{equation}

while the maximal amplitude of the potential at infinity is

(2.90) \begin{align} |\phi (\infty ,\pi )|=\lim _{r\to \infty }\max _{[r,r+2\pi /k]}|\phi (\boldsymbol{\cdot },\pi )|=2. \end{align}

The following two statements concerning phase shift and run-up can be found, for example, in Sretensky (Reference Sretensky1977) and Stoker (Reference Stoker2019) when the author discusses the known results for the case of discrete solutions $\alpha =\pi /(2m)$ , $m=1,2,\ldots$ . On the other hand, a common expectation noted in most of the classic works is that the formulas for discrete values are valid in the general case (arbitrary angle) due to continuity arguments. We formulate these statements using our notation and verify their validity numerically using our approach.

Statement 2. The reflection coefficient (2.83) for arbitrary $0\lt \alpha \leqslant \pi$ ( $\theta _0=\alpha +\pi$ ) can be calculated by the formula

(2.91) \begin{equation} R=e^{\frac {i\pi (\pi -2\alpha )}{4\alpha }}=e^{\frac {i\pi (3\pi -2\theta _0)}{4(\theta _0-\pi )}},\quad \nu =\frac {\pi (\pi -2\alpha )}{8\alpha }. \end{equation}

Statement 3. The following equation is valid for arbitrary $0\lt \alpha \leqslant \pi$ ( $\theta _0=\alpha +\pi$ ):

(2.92) \begin{equation} \frac {1}{2}\,e^{-i\nu }\phi (0,\pi )=\sqrt {\frac {\pi }{2\alpha }}=\sqrt {\frac {\pi }{2(\theta _0-\pi )}}. \end{equation}

Since Statements 2 and 3 have been theoretically proven for $\alpha =\pi /(2m)$ , we use such cases to additionally verify the accuracy of our computations.

Corollary 1. We can immediately conclude from ( 2.89 ) and ( 2.92 ) that

(2.93) \begin{equation} \tilde p_*^{(1)}(0)\, \cos \nu =\sqrt {\frac {\pi }{2\alpha }}-\sin \nu -\cos \nu = \sqrt {\frac {\pi }{2\alpha }}-\sqrt {2}\,\sin \frac {\pi ^2}{8\alpha }. \end{equation}

Note that, for all cases $\alpha =\pi /(6+8k)$ , where $k=0,1,2,\ldots$ (that is, $m=3+4k$ in terms of the discrete formulation), $ \tilde p_*^{(1)}(0)$ takes its infinite value while the multiplication of ( 2.93 ) is always finite.

The only question that remains unanswered is the behaviour of the term $T(y) \cos \nu$ in (2.88) at zero and infinity, where the results of Appendix A are instrumental. Indeed, using the residuum theorem we conclude that for all cases, except the limiting one when $\alpha =\pi$ , the following estimate is valid:

(2.94) \begin{equation} T(y)=\tilde p_*^{(1)}(0)-\frac {\pi y}{2}P_{-1}+\Gamma \left (1-\frac {\pi }{\alpha }\right )\cos \frac {\pi ^2}{2\alpha }P_*^{-}y^{\frac {\pi }{\alpha }}+ O\big (y^2\log y\big ),\quad y\to 0, \end{equation}

where $\pi /\alpha \gt 1$ . In the limiting case, $\alpha =\pi$ , this estimate changes to

(2.95) \begin{equation} T(y)=\tilde p_*^{(1)}(0)-\frac {\pi y}{2}\big (Q_*^{-}+Q_{-1}\big (1-\gamma -\log y\big )\big )+ O\big (y^2\log y\big ),\quad y\to 0, \end{equation}

where $\gamma =-\Gamma '(1)\approx 0.5772156649$ is the Euler–Mascheroni constant see (for example Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014)).

To estimate the behaviour of the function $T(y)$ at infinity we note that the poles of the function $\tilde p_*^{(1)}(s)$ with positive real parts are simple and located at points $s_1^{(l)}=2l+1$ , $l=1,2,\ldots$ . On the other hand, the function $\cos (\pi s/2)$ has zeros at those poles. This might suggest an exponential decay for the function $T(y)$ . However, there might also exist simple poles of the function $\tilde p_*^{(1)}(s)$ at points $s_2^{(l)}=1+\pi (1+2l)/(2\alpha )$ , $l=0,1,\ldots$ , and, consequently,

(2.96) \begin{equation} T(y)=\Gamma \left(1+\frac {\pi }{2\alpha }\right)\sin \frac {\pi ^2}{4\alpha } P_*^+ y^{-1-\pi /(2\alpha )}+O \big (y^{-1-3\pi /(2\alpha )}\big ),\quad y\to \infty . \end{equation}

Thus, the decay is, generally speaking, not exponential. Note that the estimate (2.96) is the worst-case scenario. It does not contradict the known classic results (see Sretensky (Reference Sretensky1977) and Stoker (Reference Stoker2019)), that, for $\alpha =\pi /(2m)$ , the corresponding coefficients $P_*^+$ are zeros. See the discussions after (A6) in Appendix A.

All constants in the asymptotic estimates (2.94)–(2.96) are defined in Appendix A. Moreover, $\tilde p_*^{(1)}(0)$ and $P_{-1}$ and $Q_{-1}$ are evaluated there explicitly, while $P_*^{-}$ , $Q_*^{-}$ and $P_*^{+}$ require computations of specific functionals acting on the solution of the integral equation.

Interestingly, taking estimate (2.94), Statements 2 and 3 and (A9) into account, the following immediately follows.

Statement 4. For arbitrarily $0\lt \alpha \lt \pi$ , there is an explicit relationship between the value of the potential $\phi (0,\pi )$ and its partial derivative over the boundary $\phi '_r(0,\pi )$ at the shoreline edge

(2.97) \begin{equation} \frac {\partial \phi }{\partial r}(0,\pi )=-2k\,e^{i\nu }\sqrt {\frac {\pi }{2\alpha }}\cot \alpha =-k \, \phi (0,\pi ) \cot \alpha . \end{equation}

To the best of our knowledge, this result is new and might be important for the analysis of wave breakage near the shore edge.

Finally, when $\tilde \phi (s,\pi )$ has been found, we can also reconstruct the entire solution inside the water from (2.20)

(2.98) \begin{align} \phi (r,\theta )& =\frac {1}{2\pi i}\int _{-i\infty +\eta }^{i\infty +\eta } \frac {\Gamma (s)}{k^{s}}\frac {\cos s(\theta -\theta _0)}{\cos s(\pi -\theta _0)} \nonumber \\[3pt] & \quad \times\left ((1+R)\big (1+\tilde p^{(1)}_*(s)\big )\cos \frac {\pi s}{2}-(1-R)i\sin \frac {\pi s}{2}+i(1-R)\right ) r^{-s}\text{d}s. \end{align}

Here, $\eta$ can be taken equal to zero, while all functions decay exponentially along the imaginary axis. Moreover, a solution on the inclined boundary ( $\theta =\theta _0$ ) can be written in the form

(2.99) \begin{equation} \phi (r,\theta _0)=\frac {1}{2\pi i}\int _{-i\infty +\eta }^{i\infty +\eta } \frac {\tilde \phi (s,\pi )}{\cos s\alpha }r^{-s}\text{d}s=\frac {1}{\alpha } \int _0^\infty \frac {(rt)^{\frac {\pi }{2\alpha }}\, \phi (t,\pi )} {t^{\frac {\pi }{\alpha }}+r^{\frac {\pi }{\alpha }}}\frac {\text{d}t}{t}, \end{equation}

since the solution on the surface of the free water has already been computed in (2.88). Note that $\phi (r,\theta _0)\to \phi (r,\pi )$ as $\theta _0\to \pi$ . Another trivial fact is $ \phi (0,\theta _0)= \phi (0,\pi ).$ For numerical computations of (2.99), the following identity can be useful:

(2.100) \begin{align} \int _0^\infty \frac {(rt)^{\frac {\pi }{2\alpha }}\, f(t)} {t^{\frac {\pi }{\alpha }}+r^{\frac {\pi }{\alpha }}}\frac {\text{d}t}{t}=\frac {\alpha f(0)}{1+r^{\frac {\pi }{2\alpha }}}+ \int _0^\infty \frac {(rt)^{\frac {\pi }{2\alpha }}\, } {t^{\frac {\pi }{\alpha }}+r^{\frac {\pi }{\alpha }}}\left (f(t)-\frac {f(0)}{1+t^{\pi /\alpha }}\right )\frac {\text{d}t}{t}. \end{align}

3.1. Numerical scheme to solve the integral equation (2.77)

In our calculations, we replace the integral defined on the half-axis by the one on a finite domain $(\epsilon ,1/\epsilon )$ . We use a logarithmic mesh and Simpson’s composite quadrature rule (exploiting the smoothness of the solution). Near the point $r=1$ , we take into account that the integral is understood in terms of Cauchy singular integration. The asymptotic behaviour (2.80) discussed above can be used to control the accuracy of the computation.

As a result, we compute a solution of the respective integral equation as a function of the parameter $d_1$ . The solution, generally speaking, has a singularity of this kind

(3.1) \begin{equation} p_*^{(1)}(r,d_1)\sim \frac {\varpi (d_1)}{r-1}=\frac {\varpi _1-d_1\varpi _d}{r-1} ,\quad r\to 1, \end{equation}

that can be easily traced numerically. We optimise the solution with respect to the parameter $d_1$ to ensure its smoothness at point $r=1$ by requiring the condition $\varpi (d_1)=0$ . Since $\varpi (d_1)$ is a linear function of the parameter $d_1$ , one can use any method to compute its unique value ( $d_1=\varpi _1/\varpi _d$ ).

3.2. Comparison of the numerical solution with the exact one for discrete case: $\alpha =\pi /(2m)$

Note that (2.77) can be solved analytically in the discrete case.

The case $\theta _0=\pi +\pi /2$ . We observe that $n_0(r,\theta _0)\equiv 0$ and thus the integral operator $\mathcal N$ vanishes. As a result, the unique solution of the integral equation is $q_*^{(1)}=0$ and $d_1=0$ . Due to (2.82), we conclude that $R=1$ , and the entire solution is represented by the function $U_0(r)$ describing waves at infinity, as expected.

The case $\theta _0=\pi +\pi /4$ . Now we can check that $n_0(r,\theta _0)= -n_d(r,\theta _0)$ . Thus, the only solution from the required class is $d_1=-1$ , $p^{(1)}_*(r)=0$ and $R=i$ (compare with (2.82)) or the phase shift, $\nu =\pi /4$ , (2.83). This result naturally coincides with the well-known one from Lewy (Reference Lewy1946), Stoker (Reference Stoker1947) and Sretensky (Reference Sretensky1977).

The case $\theta _0=\pi +\pi /(2m)$ . For those cases $n_0(r,\theta _0)$ and $n_d(r,\theta _0)$ are no longer linearly dependent. However, if one specifies the representation (2.60) even further (using the arguments $\pi s/(2m)$ , $m=2,3,4 \ldots M$ in the trigonometric functions), the respective function $n_d$ will change and the solution of the equation can still be computed in a close form, while the solution will be the same cumbersome as in the classic case (compare with Lewy (Reference Lewy1946)). Those solutions will decay exponentially at infinity and the respective phase shift is equal to $\nu =\pi (m-1)/4$ , as expected (Lewy Reference Lewy1946; Stoker Reference Stoker1947). We do not go in this direction; simply compute the integral equation for some of the discrete angles $\theta _0=\pi +\pi /(2m)$ .

We start with the verification of our integral equation (2.77) (and the entire solution) comparing it with the discrete cases $\alpha =\pi /(2m)$ , $m=1,2,\ldots$ . Note that the solution of the integral equation for the cases $m=1,2$ should be, as discussed above, identically zero. In fact, the respective curves in figure 2(a) lie below the threshold of $10^{-15}$ , which is the precision of the MATLAB calculations. Note that the spikes observed in the curves plotted in log–log scale correspond to the zeros of the solution of the integral equation, that is, to the points where these real functions change sign.

The remaining six curves behave as expected in the estimate (2.80) with the exponents $\beta _0=1$ and $\beta _\infty =3$ . Those facts are consistent with the analysis provided in Appendix A (compare with discussion around (A5)).

Figure 2.

(a) Absolute value of the (real) solutions $p_*^{(1)}(t)$ of the integral equation (2.77) for different angles $\alpha =\theta _0-\pi =\pi /(2m)$ , $m=1,2,\ldots ,8$ . The graphs are presented in log–log scale thus the spikes inside the interval correspond to zero values of the solutions. (b) Shows the respective inversion of the second term in the representation (2.88). Note that the curves for $m=1,2$ are identically zero, while those for $m=3,7$ are close to zero due to the extremely small value of the multiplier, $\cos \nu$ .

Note that the values of the solution of the integral equation, $p_*^{(1)}(t)$ , for $m=3,7$ are a few orders of magnitude higher than for the other cases (compare with figure 2(a)). In fact, if one reformulates the problem with new unknown functions $\tilde p_*^{(j)}(s)\cos \nu$ instead of $\tilde p_*^{(j)}(s)$ , it would work for all cases.

To make this point exposed, in figure 2(b) we presented the inverted terms in the representation (2.88). It represents the main part of the localised part of the solution. Evanescent waves are drawn in figure 4(a), while the entire water surface profile is given in figure 4(b).

Having the function $p_*^{(1)}(t)$ at hand, we can compute the value $\phi (0,\pi )$ from (2.88) and thus verify the identity (2.92), which has been proved valid for any discrete solution. Note that the values of the reflection coefficients $R$ are computed from the formula (2.82) when solving the integral equation (see the discussion around (3.1). The respective results are presented in figures 3(a) and 3(b).

We observe that as the angle decreases, the computational challenges increase. This has been predicted (see discussions around (2.78)). It is clear that, by increasing the number of mesh points, we could improve the accuracy if needed. However, for extremely small angles $\alpha$ , the shallow-water approximation should be used naturally instead. Moreover, the sensitivity of the reflection coefficient becomes an issue and should be questioned against the physical phenomena behind the potential flow.

Figure 3.

Indirect verification of the computational accuracy of the integral equation (2.77) for different angles $\alpha =\theta _0-\pi =\pi/(2m)$ , $m=1,2,\ldots ,8$ by comparing the numerical and exact analytical values of the reflection coefficient and the value of the potential at zero point along the water surface ( $z=0$ ).

Figure 3. Long description

Panel A: A scatter plot shows the relationship between alpha over pi and the absolute value of phi at zero and pi. The x-axis is labeled alpha over pi equals theta zero over pi minus one equals one over two m, with values ranging from 0 to 0.5. The y-axis is labeled absolute value of phi at zero and pi, with values ranging from 2.0 to 6.0. The plot includes data points in various colors and an inset graph showing relative error versus alpha over pi, with values ranging from 0 to 0.6 on the x-axis and 0 to 1 times 10 to the power of negative 3 on the y-axis. Panel B: A line graph shows the e^(-iν)φ(r,π) as a function of negative k times r. The x-axis is labeled negative k times r, with values ranging from -10 to 0. The y-axis is labeled e^(-iν)φ(r,π), with values ranging from -3 to 7. The graph includes multiple lines representing different values of m, ranging from 1 to 8, with a legend indicating the corresponding colors and line styles.

Figure 4.

(a) Evanescent waves along the water surface for different angles $\alpha =\theta _0-\pi =\pi/(2m)$ , $m=1,2,\ldots ,8$ . (b) Shows the full water surface. In case $m=1,2$ , the solution is given analytically (respective curves are given by dashed and dotted lines.

In figure 4(a) we present the localised part (evanescent waves), while in figure 4(b) we show the entire profile of water surface ( $z=0$ ) for those values of the angle $\alpha =\pi /(2m)$ , $m=1,2,\ldots$ . In conclusion of this subsection, we confirm that the new method produces the same results as known from Lewy (Reference Lewy1946), Stoker (Reference Stoker1947) and Sretensky (Reference Sretensky1977). Note that these cases correspond only to the reflection coefficients $R=\pm 1,\pm i$ .

3.3. Verification of the Statements 2–3 for arbitrary angles $0\lt \alpha \leqslant \pi /2$

In this section, we compute the solution in the same region $0\lt \alpha \lt \pi /2$ . In our computations, we can take an arbitrary number. On the other hand, any angle can always be represented with an arbitrary accuracy in the form $\alpha =\pi p/(2m)$ and $p\lt 2m$ where the theory is available, but the respective analytical solution takes too complicated form Lewy (Reference Lewy1946), Stoker (Reference Stoker1947) and Sretensky (Reference Sretensky1977). By the continuity argument, the global parameters of the solution as the phase shift (reflection coefficient) and the run-up ratio should be the same for angles $\alpha \lt \pi /2$ . We verify below the Statements 2 and 3 for $\alpha =\pi /(2l+1)$ , $l=1,2,\ldots$ .

Figure 5.

(a) Absolute value of the solutions $p_*^{(1)}(t)$ of the integral equation (2.77) for different angles $\alpha =\theta _0-\pi =\pi /(2l+1)$ , $l=1,2,\ldots ,6$ . The graphs are presented in log–log scale thus the spikes inside the interval correspond to zero values of the solutions. Note that the value of the exponent $\beta _0=1$ for all cases. Differently for the exponent defining decay at infinity. Here $\beta _\infty =3$ for all cases except $l=1$ where $\beta _\infty =1+\pi /(2\alpha )=5/2\lt 3$ (see for details in Appendix A). An arrow depicts this respective curve. (b) Shows the respective inversion of the second term in the representation (2.88).

In figure 5(a) we present the absolute value of the numerical solutions of the integral equation (2.77) on the log–log scale. The only essential change compared with the results of the previous subsection (figure 2 a) is the fact that for $\alpha =\pi /3$ , the behaviour at infinity from the estimate (2.80) changes from $\beta _\infty =3$ to $\beta _\infty =1+\pi /(2\alpha )=3/2\lt 3$ . Such an effect has been predicted in the more detailed analysis of the solution given in Appendix A. It was not previously observed for the angles $\alpha =\pi /(2m)$ , and the values $\alpha =\pi /4$ and $\pi /2$ led to a trivial solution of the integral equation. We also note that in the cases under consideration all solutions have maximal values of the order of $10^0$ instead of $10^3$ as previously for those cases $m=3,7$ where $d_1$ were large.

In figure 5(b), as previously, we present the second term in the representation (2.88) for different values of the angle $\alpha$ . Now the multiplier $\cos \nu$ plays no role in the calculations as it always takes the same absolute value 1/ $\sqrt {2}$ .

Figure 6.

Verification of the Statements 2 and 3 for different angles $\alpha =\theta _0-\pi =\pi/(2l+1)$ , $l=1,2,\ldots ,6$ by comparing the numerical and the predicted analytical results.

Figure 6. Long description

Panel A: A scatter plot shows the magnitude of φ(r,π) on the vertical axis against α/π on the horizontal axis. The data points are marked with different symbols and colors, indicating various series. An inset within this panel displays a smaller scatter plot with a different scale for the vertical axis, showing relative error against the same horizontal axis ratio. Panel B: Another scatter plot compares the reflection coefficient, with real and imaginary parts plotted against the same horizontal axis ratio. The legend indicates different symbols and colors representing real and imaginary parts of both numerical and analytical results. The data points are clustered at specific values, and annotations indicate where all values coincide.

Figure 7.

(a) Evanescent waves along the water surface for different angles $\alpha =\theta _0-\pi =\pi/(2l+1)$ , $l=1,2,\ldots ,6$ . (b) Shows the full water surface.

Figure 7. Long description

Two line graphs illustrate wave behavior for different angles. Panel A shows the graph of T(kr)cos(ν) + e^(-kr)sin(ν) with the horizontal axis labeled as -kr and the vertical axis ranging from -2 to 4. The legend indicates different lines for l values ranging from 1 to 6, each showing distinct trends. Panel B displays the graph of e^(-iν)φ(r,π) with the horizontal axis labeled as -kr and the vertical axis ranging from -3 to 6. The legend similarly indicates different lines for l values ranging from 1 to 6, each exhibiting unique patterns.

In figures 6(a) and 6(b) we verify the formulas of Statements 2 and 3. The accuracy is on the same level (or better) as in the previous subsection. Finally, in figures 7(a) and 7(b), we present the evanescent waves and the full profile of the potential along the water surface ( $z=0$ ).

Thus, we can confidently conclude that Statements 2 and 3 are valid for the angles $0\lt \alpha \leqslant \pi /2$ . The deviations of the numerically computed values from the analytical formulas are due to the numerical error. The latter increases for smaller angles of $\alpha$ as we have already discussed in the previous subsection.

In the next subsection, we will analyse the formulas (2.92) and (2.91) from Statements 2 and 3 that represent the phase shift and the run-up for the values of $\alpha$ from the interval $\pi /2\lt \alpha \leqslant \pi$ . We use the same numerical scheme to solve the integral equation (2.77).

3.4. Verification of the Statements 2–3 for arbitrary angles $\pi /2\lt \alpha \leqslant \pi$

In figure 8(a), we show the respective solutions of the integral equations for a few values of ( $\alpha /\pi =0.51,0.6,0.7,0.8,0.9,1.0$ ). As seen, now the exponent, $\beta _\infty =\beta _\infty (\alpha )$ in the asymptotics at infinity (2.80) depends on the angle $\alpha$ , while $\beta _0=1$ remains. Only in the limiting case ( $\alpha =\pi$ ), the leading term in the asymptotics (2.80) $_1$ changes to $t\log t$ as predicted in Appendix A. This fact is hardly visible, so we indicate it by an arrow.

Figure 8.

(a) Represents the absolute value of the solutions $p_*^{(1)}(t)$ of the integral equation (2.77) for different angles $\alpha =\theta _0-\pi$ as indicated in the legend. The graphs are presented in log–log scale thus the spikes inside the interval correspond to zero values of the solutions. The value of the exponent $\beta _0=1$ for all cases except for $\alpha =\pi$ when the asymptotic is different (it contains also a logarithmic term). As for the other exponent, $\beta _\infty$ , it now depends on the angle $\beta _\infty =1+\pi /(2\alpha )\in (3/2,2)\lt 3$ (see for details in Appendix A). (b) Shows the respective inversion of the second term in the representation (2.88).

The second term in the solution representation (2.88) now changes its behaviour. It no longer indicates the oscillatory character in contrast to the case $0\lt \alpha \lt \pi /2$ . Also, the values are an order of magnitude smaller than in the case $0\lt \alpha \lt \pi /2$ while the decay is different (no longer exponential).

Figure 9.

Verification of the formulas for the run-up and the reflection coefficient as shown in the legend from Statements 2 and 3 for different angles $\alpha \gt \pi /2$ . For the same numerical settings, the computational accuracy is higher than for smaller angles $\alpha \lt \pi /2$ .

Figure 9. Long description

Panel A: A line graph shows computations and a theoretical function e^(-ivφ) plotted against α/π. The x-axis is labeled α/π ranging from 0.5 to 1.0, and the y-axis ranges from 1.4 to 2.0. Computational data points are marked with circles, and the theoretical function is represented by a solid line. An inset graph titled ‘Relative error’ shows a scatter plot of relative error against α/π, with the x-axis ranging from 0.5 to 1.0 and the y-axis ranging from 10^-6 to 10^-4. Panel B: A line graph compares the reflection coefficient R for different angles. The x-axis is labeled α/π = θ0/π-1 ranging from 0.5 to 1.0, and the y-axis ranges from -0.8 to 1.0. The graph includes lines for sin 2v, cos 2v, Im R, and Re R, with corresponding markers and colors indicated in the legend. An annotation points to a specific data point labeled ‘Dock problem’.

Figure 10.

(a) Localised part of the solutions along the water surface for different angles $\alpha$ as indicated in the legend, while (b) shows the full profile of potential at $z=0$ .

In figures 9(a) and 9(b), we verify the formulas for the reflection coefficient and the run-up from Statements 2 and 3 for different angles $\alpha$ as shown in the legend. For the same numerical settings, the computational accuracy is higher than for smaller angles $\alpha \in (0,\pi /2)$ . Moreover, the variation of the reflection coefficient is much slower now (no oscillations at all). In the limiting case ( $\alpha =\pi$ ) the values coincide with those of the dock solution (Chakrabarti et al. 2005).

To conclude analysis of this case, we present in figures 10(a) and 10(b), the entire localised part of the solution near the full water profile along the water surface, respectively. As mentioned above, the localised part is much less pronounced for $\pi /2\lt \alpha \leqslant \pi$ than in the case of smaller angles, that is, $0\lt \alpha \lt \pi /2$ .

To summarise the verification from Statements 2 and 3 in the entire range of the beach slope $0\lt \alpha \leqslant \pi$ , in figure 11, we present analytical and numerical results for the reflection coefficient and the run-up at the beach edge.

Figure 11.

Summary on the validity of the formulas from Statements 2 and 3 in the entire domain of the angle $0\lt \alpha \leqslant \pi$ . All available results have been used for the verification. Accuracy of the computations for the fixed set of technical parameters when computing the solution $p^{(1)}_*(t)$ of the integral equation and evaluating the function  $T(y)$ .

3.5. Computation of pressure on the seabed for various beach angles $0\lt \alpha \leqslant \pi$

In this subsection, we analyse the impact of water waves on the seabed, with the pressure computed as $ p(r,t)=p_0+\rho g r\sin \alpha +\mathrm{Re} \{i\rho \omega \phi (r,\theta _0)e^{-i\omega t}\}$ or

(3.2) \begin{align} p(r,t)=p_0+\rho g r\sin \alpha +\rho \omega \cos \left (\frac {\pi }{2}+\nu -\omega t\right )e^{-i\nu }\phi (r,\theta _0). \end{align}

In figure 12(b), we present the distribution of the normalised potential $e^{-i\nu }\phi (r,\theta _0)$ along the water bed ( $\theta =\theta _0=\alpha +\pi$ ) for various angles $\alpha$ smaller than $\pi /2$ (figure 12 a) and larger than $\pi /2$ (in figure 12(b)). It preserves the same features in terms of monotonicity as the free water surface.

Figure 12.

The figure presents the normalised potential along the water bed ( $\theta =\theta _0=\alpha +\pi$ ) for various angles $\alpha$ smaller than $\pi /2$ (a) and larger than $\pi /2$ ((b)).

Figure 12. Long description

Panel A: A line graph shows the normalized potential along the water bed for various angles smaller than π/2. The horizontal axis is labeled ‘kr’ and the vertical axis is labeled ‘e^(-iv)φ(r, θ0), θ0 = α + π’. The legend indicates different angles (α) with corresponding colors: α = π/13 (blue), α = π/11 (orange), α = π/9 (yellow), α = π/7 (purple), α = π/5 (green), α = π/3 (light blue), and α = π/2 (black). The lines show oscillatory behavior with varying amplitudes and frequencies. Panel B: Another line graph shows the normalized potential along the water bed for various angles larger than π/2. The horizontal axis is labeled ‘kr’ and the vertical axis is labeled ‘e^(-iv)φ(r, θ0), θ0 = α + π’. The legend indicates different angles (α) with corresponding colors: α = 0.52π (blue), α = 0.6π (orange), α = 0.7π (yellow), α = 0.8π (purple), α = 0.9π (green), α = 1.0π (light blue), and α = π/2 (black). The lines show a decreasing trend with varying slopes.

3.6. Numerical validation using the boundary element method

Although we have strong evidence of the accuracy of our results, we also validated them using the BEM. This allowed us to check the reflection coefficient and the free surface displacement. The BEM is a numerical method that is ideally suited to solving the Laplace equation. For our problem, we use a modification of the method outlined in Wang & Meylan (Reference Wang and Meylan2002). Note that the BEM code can be downloaded from WikiWaves (2010). Figure 13 shows the finite triangular region used to approximate the solution with the boundary conditions. All details can be found in Wang & Meylan (Reference Wang and Meylan2002). Note that the left-hand boundary condition is not pointwise because it involves a linear operator $\mathbf{Q}$ . For our calculations, we assume $k=1$ and the depth (of the computational domain) is $h=4$ .

Figure 13.

Geometric setting used for the BEM computations. The boundaries are the dashed lines with the normal derivative given by the equations shown.

The BEM that we use has linear elements and the convergence will be linear with the size of the panel. We consider two panel sizes, approximately 0.025 and 0.0125 (note that these are by necessity approximate since the lengths of each sides will not necessarily be divisible by the panel size). Figure 14 shows the exact and BEM solution for the reflection coefficient, $R$ , as a function of $\alpha$ and $\nu$ . As expected, the error increases as the angle decreases and behaves approximately as $1/\alpha$ for small values of $\alpha$ . Note that the same panel size is used here, so that as the angle increases, the size of the BEM matrices increases significantly. Figure 15 shows the solution for the two panel lengths. On the left, we show the error of the reflection coefficient and, on the right, the solution for the displacement at the end. Note that the end displacement is calculated using a linear extrapolation of the second and third last elements because there are large errors for the final element. Linear convergence with panel size is apparent, as is the linear growth in error with respect to $\nu$ (we note that this linear growth with $\nu$ does not follow from the linearity of the elements).

Finally, we note that the BEM is a numerical method with great flexibility, for example, it can be extended to non-uniform beach profiles, obstacles etc. For more complicated problems, for example, where there is a singularity at the end point, and for which no analytic results exist, the BEM would provide a useful validation.

Figure 14.

Real and imaginary part of $R$ as a function of $\alpha /\pi$ (a) and $\nu$ (b). The solid line is the exact solution and the dots are the solution by the BEM, with panel size approximately $0.025$ .

Figure 14. Long description

Panel A: A line graph shows the real and imaginary parts of a function as a function of alpha over pi. The x-axis is labeled alpha over pi, ranging from 0 to 1. The y-axis is labeled Re (R) and Im (R), ranging from -1.0 to 1.0. The graph includes blue stars representing the real part of R, red stars representing the imaginary part of R, a blue line representing cos 2v, and a red line representing sin 2v. The real and imaginary parts of R oscillate and intersect with the cosine and sine functions. Panel B: A line graph shows the real and imaginary parts of a function as a function of nu. The x-axis is labeled nu, ranging from 0 to 20. The y-axis is labeled Re (R) and Im (R), ranging from -1.0 to 1.0. The graph includes blue stars representing the real part of R and red stars representing the imaginary part of R. The real and imaginary parts of R exhibit oscillatory behavior.

Figure 15.

The absolute value of the difference between the reflection coefficient, $|R|$ , calculated by the BEM and the exact solution (a) and the absolute value of the potential at the end, $|\phi (0,\pi )|$ , calculated using the BEM and the exact formula (b). Note that the step sizes are approximately $0.025$ (red) and $0.0125$ (green).

4.1. Numerical calculation

Solving (4.2) numerically can be accomplished straightforwardly by converting to matrix multiplication as outlined in Wilks et al. (Reference Wilks, Meylan, Montiel and Wakes2024, Reference Wilks, Meylan, Montiel and Wakes2025). We divide $r$ and $k$ into a set of even points $r_i$ and $k_{\!j}$ with spacing $\Delta r$ and $\Delta k$ respectively. For our numerical calculations, we choose 1000 points for $x$ from $-20$ to $0$ and for $k$ we choose 5000 points between $0$ and $5$ . We define the matrix

(4.4) \begin{equation} \mathbf{P} = [P_{ij}], P_{ij} = \phi(r_i, \pi, k_j). \end{equation}

Then

(4.5) \begin{equation} \mathrm{vector}[\zeta (r_i,t+t_0)] \approx \frac {\Delta r \Delta k}{2\pi } \mathbf{P}\times \mbox{diag}[\cos(\omega_jt)]\times \mathbf{P}^{\star } \times \mathrm{vector}[\zeta (r_i,t_0)], \end{equation}

where vector means to convert the sequence of points to a column vector and diag denotes the same conversion to a diagonal matrix. The star denotes the complex conjugate transpose. We note that we only need to calculate the matrix $\mathbf{P}$ once and can then solve for any initial condition and at any time. Except for the fact that we do not consider the initial value of the time derivative (which is only to keep the explanation simple), we have effectively calculated the operator $\exp (-\mathrm{i} \mathcal{P} t )$ and expressed it as matrix multiplication.

We also note that we could choose the number of points in $r$ and $k$ to be the same, say, $N$ and the spacings $\Delta r$ and $\Delta k$ in such a way that

(4.6) \begin{align} \frac {\Delta r \Delta k}{2\pi } = \frac {1}{N}. \end{align}

This would in some sense be the optimal choice, and we note that this is what is done when converting the Fourier transform to the fast Fourier (or discrete) Fourier transform in the standard implementations. Further details of this can be found in Wilks et al. (Reference Wilks, Meylan, Montiel and Wakes2024, Reference Wilks, Meylan, Montiel and Wakes2025).

We present some time-domain results with a solution with an initial condition

(4.7) \begin{equation} \zeta_0 = 2e^{-(x+3)^2}, \quad x\le0,\quad (r=-x,\quad \theta=\pi). \end{equation}

Figure 16.

The evolution of the fluid surface for the initial condition given by (4.7) with no initial velocity. The slopes are $\alpha = 0.089\pi , 0.186\pi$ and $0.380\pi$ . The non-dimensional times (gravity is assumed unity) are shown. The animations can be found in the supplementary material (Movies 1–3), available at https://doi.org/10.1017/jfm.2026.12121.

Figure 16. Long description

The image contains eighteen line graphs arranged in a 6x4 grid. Each graph shows the evolution of the fluid surface over time for different slope angles and non-dimensional times. The x-axis represents the spatial coordinate (x) ranging from -20 to 0, and the y-axis represents the fluid surface elevation (ζ) ranging from -3 to 3. Each row corresponds to a different non-dimensional time (t = 0, t = 10, t = 20, t = 30), and each column corresponds to a different slope angle (α = 0.089π, α = 0.186π, α = 0.380π). The graphs illustrate how the fluid surface evolves over time for these specific conditions. Panel A: The first row shows the initial condition at t = 0 for three different slope angles. Panel B: The second row shows the fluid surface at t = 10 for the same slope angles. Panel C: The third row shows the fluid surface at t = 20. Panel D: The fourth row shows the fluid surface at t = 30. Each graph captures the changes in the fluid surface elevation over the specified spatial range and time intervals.

Figure 16 shows the evolution of the fluid surface for the initial condition given by (4.7) with no initial velocity. The slopes are $\alpha = 0.089\pi , 0.186\pi$ and $0.380\pi$ . The non-dimensional times (gravity is assumed unity) are shown. The animation can be found in the supplementary material as Movies 1–3. Figure 16 shows the evolution of the slopes $\alpha = 0.089\pi$ , $0.186\pi$ and $0.380\pi$ , respectively. Figure 17 is the same as figure 16 except for $\alpha = 0.477\pi , 0.670\pi$ and $0.864\pi$ . The animation for this figure can be found in the supplementary material as Movies 4–6. Note that the solution for these and subsequent figures is given for non-dimensional time (assuming gravity is unity). The pulse initially travels in both directions and the right-travelling wave leads to a high amplitude on the beach before being reflected. As the angle of the beach increases the run-up decreases, and also the figures become more and more similar, reflecting the fact that the variation in the frequency-domain solutions becomes less as the angle increases.

Figure 17.

As in figure 16 except $\alpha = 0.477\pi , 0.670\pi$ and $0.864\pi$ . The animations can be found in the supplementary material as Movies 4–6.

4.2. Solution for an incident wave packet

The formulation given expresses the solution as a weighted integral of the frequency-domain solutions, where the weighting is calculated using an inner product. A simpler approach is to impose the weighting and express the solution as

(4.8) \begin{equation} \zeta (r,t) =\mathrm{Re}\left \{\int _{0}^{\infty } \hat {f}(k) \phi (r,\pi ,k) \exp (-\mathrm{i} \omega t) \mathrm{d}k\right \}\!, \end{equation}

where we specify $\hat {f}(k)$ . For the numerical examples, we choose

(4.9) \begin{equation} \hat {f}(k) = \sqrt {\frac {5}{\pi }}\exp \big (-5(k-2)^2)\big ). \end{equation}

Note that this calculation can be transformed to a matrix multiplication, as outlined in the previous section. This is not an initial value problem. The wave packet is incident from negative infinity and $t=0$ is set to be the time the centre of the pulse would reach $r=0$ if the depth remained infinite. Therefore, we can consider negative time.

Figure 18 shows the evolution of the wave packet given by (4.9) for the non-dimensional times (gravity is unity) shown with slopes $\alpha = 0.089\pi , 0.186\pi$ and $0.380\pi$ . The animation can be found in the supplementary material as Movies 7–9. Figure 19 is the same as figure 18 except for $\alpha = 0.477\pi , 0.670\pi$ and $0.864\pi$ . The animation for this figure can be found in the supplementary material as Movies 10–12. As shown in figures 16 and 17 the variation in solution is less as the angle increases. The amplification of the pulse at the beach end is clearly apparent for small angle $\alpha$ .

Figure 18.

The evolution of the wave packet given by (4.9) for the non-dimensional times (gravity is unity) shown with slopes $\alpha = 0.089\pi , 0.186\pi$ and $0.380\pi$ . The animations can be found in the supplementary material as Movies 7–9.

Figure 18. Long description

The image contains twelve line graphs arranged in a 4x3 grid, each showing the evolution of a wave packet over time for different slope angles and non-dimensional times. Each graph has the horizontal axis labeled ‘x’ ranging from -20 to 0 and the vertical axis labeled ‘ξ’ ranging from -3 to 3. The graphs are organized by rows and columns, with each row corresponding to a specific non-dimensional time (t = -20, t = -10, t = 0, t = 10, t = 20, t = 30) and each column corresponding to a specific slope angle (α = 0.089π, α = 0.186π, α = 0.380π). Each graph shows a blue line representing the wave packet’s evolution. The wave patterns vary across different times and slope angles, illustrating how the wave packet changes over time and with different slope conditions.

Figure 19.

As in figure 18 except $\alpha = 0.477\pi , 0.670\pi$ and $0.864\pi$ . The animations can be found in the supplementary material as Movies 10–12.

Finally, we compute the wave height on the shoreline $\zeta (t,0)$ divided by $\tan ^2(\alpha )$ for the initial pulse (figure 20 a) and the incident pulse (figure 20 b). The values for all parameters are the same as in figures 16 and 18. The normalisation of this kind has been chosen to make both results equally visible. It was suggested by Miche’s wave breakage condition (Miche Reference Miche1944) that states wave breaking occurs when the slope of the seabed, $S$ , becomes smaller than the gradient of the wave-induced surface elevation, $ \max _t\zeta '_x(t,0) \gt S=\tan \alpha ,$ where $\zeta (t,x)$ is the elevation of the water surface, computed from (4.2) for the initial pulse and from (4.8) for the incident packet. For monochromatic waves, Miche’s criterion in our notation (compare with (2.97)) reads: $ \max _t \zeta (t,0)\gt k^{-1}\tan ^2\alpha .$ The mixture of frequencies in these two cases is very different, but the scaled amplitude on the beach in both cases is similar. Since the initial pulse breaks into two waves which travel in opposite directions, in some sense both waves have non-dimensional amplitude unity near the beach.

Figure 20.

The height of the wave at the shoreline $\zeta (t,0)$ divided by $\tan ^2(\alpha )$ for the initial pulse (a) and the incident pulse (b). The values for all parameters are the same as in figures 16 and 18.