2.1. Equations of motion and useful relations

We consider the steady two-dimensional motion of an incompressible inviscid fluid with constant vorticity $\omega$. The fluid is bounded above and below by an impermeable free surface and solid horizontal bed, respectively. Our focus lies on travelling waves of permanent form that propagate with a constant phase speed $c$ and wavenumber $k$ ($k=0$ for solitary and more general aperiodic waves). We adopt a Galilean frame of reference moving with the wave, thus ensuring that the velocity field appears independent of time for the observer. Consequently, we can express the fluid domain, denoted as $\varOmega$, as the set of points $(x,y)$ (Cartesian coordinates) satisfying $x \in \mathbb {R}$ and $-d \leqslant y \leqslant \eta (x)$, where $\eta (x)$ represents the surface elevation from rest and $d$ is the mean water depth. Thus, the mean water level is located at $y=0$, such that

(2.1)\begin{equation} \langle \eta \rangle \stackrel{\text{def}}= \frac{k}{2{\rm \pi}} \int_{-{\rm \pi}/k}^{{\rm \pi}/k} \eta(x)\,\textrm{d}\,x = 0, \end{equation}

where $\langle \boldsymbol {\cdot }\rangle$ denotes the Eulerian averaging operator (c.f. figure 1). In the case where the surface overturns, $\eta$ is not a graph so the mapping $x\mapsto \eta (x)$ is no longer valid. We introduce in subsequent sections a parameterisation of the surface to allow for the computation and reconstruction of surface wave profiles in this context.

Figure 1. Definition sketch in the referential moving with the wave.

In this setting the velocity field $\boldsymbol {u}=(u,v)$ and pressure $p$ (divided by the density and relative to the reference value at the surface) are governed by the stationary Euler equations

(2.2a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0,\quad \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} p + \boldsymbol{g} = 0, \end{equation}

where $\boldsymbol {g}=(0,g)$, the acceleration due to gravity $g>0$ acting downwards. Equations of motion (2.2a,b) are supplemented with kinematic and dynamic conditions at the upper and lower boundaries

(2.3a)$$\begin{gather} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta = 0 \quad \textrm{at}\ y=\eta(x), \end{gather}$$

(2.3b)$$\begin{gather}p = 0 \quad \textrm{at}\ y=\eta(x), \end{gather}$$

(2.3c)$$\begin{gather}\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = 0 \quad \textrm{at}\ y={-}d, \end{gather}$$

where $\boldsymbol {n}$ is the outward normal vector (see figure 1 for a sketch of this configuration).

Since our physical system is two-dimensional per se, we introduce a scalar streamfunction $\psi$ such that $u=\psi _y$ and $v=-\psi _x$, so (2.2a) is satisfied identically and $\omega =-\psi _{xx} -\psi _{yy}$ is assumed constant. Thus, the Euler equations can be integrated into the Bernoulli equation

(2.4)\begin{equation} 2 ( p + gy ) + u^2 + v^2 = {B}_{s} - 2\omega(\psi-{\psi}_{s})\stackrel{\text{def}}= B(\psi) \end{equation}

for a constant ${B}_{s}$ (Clamond et al. Reference Clamond, Labarbe and Henry2023). Alternatively, we can define a Bernoulli constant at the bottom ${B}_{b}\stackrel {\small {\rm def}}={B}_{s} - 2\omega ({\psi }_{b}-{\psi }_{s})$. In (2.4), as in the rest of the paper, subscripts ‘s’ and ‘b’ denote that the fields are evaluated, respectively, at the surface and at the bottom. The free surface and the seabed being both streamlines, ${\psi }_{s}$ and ${\psi }_{b}$ are constant in this problem.

Because here ${p}_{s}=0$ (constant atmospheric pressure set to zero without loss of generality), we have the following relations relating some average bottom quantities and parameters (Clamond et al. Reference Clamond, Labarbe and Henry2023):

(2.5a–c)\begin{equation} \langle{p}_{b}\rangle=gd,\quad \langle{u}_{b}^{2}\rangle={B}_{b},\quad \langle{u}_{b}-(1+\eta_x^{2}){u}_{s}\rangle=\omega d. \end{equation}

Although we decide to set the reference frame as moving with the wave celerity, it is still useful to consider Stokes’ first and second definitions of phase speed,

(2.6)$$\begin{gather} c_1 \stackrel{\text{def}}={-} \langle {u}_{b} \rangle ={-}\omega d - \langle (1+\eta_x^{2}) {u}_{s} \rangle, \end{gather}$$

(2.7)$$\begin{gather}c_2 \stackrel{\text{def}}={-} \left\langle \frac{1}{d} \int_{{-}d}^\eta u \,\textrm{d} y \right\rangle = \frac{{\psi}_{b} - {\psi}_{s}}{d} ={-} \frac{\omega d}{2} - \frac{\omega\langle\eta^2\rangle}{2d} - \frac{\langle (1+\eta_x^{2}) h {u}_{s} \rangle}{d}, \end{gather}$$

where $h\stackrel {\small {\rm def}}=\eta (x)+d$ is the local water depth.

2.2. Holomorphic functions

The vorticity being constant, potential theory still holds when using a Helmholtz representation to subtract the contribution of the linear shear current (Clamond et al. Reference Clamond, Labarbe and Henry2023). Therefore, it is convenient to introduce the complex relative velocity

(2.8)\begin{equation} W(z) \stackrel{\text{def}}= U(x,y) - \textrm{i} V(x,y) = (y+d)\omega + u(x,y) - \textrm{i} v(x,y), \end{equation}

that is holomorphic in $\varOmega$ for the complex coordinate $z\stackrel {\small {\rm def}}= x+\mathrm {i} y$. (Obviously, the complex velocity $w\stackrel {\small {\rm def}}= u-\mathrm {i} v$ is not holomorphic if $\omega \neq 0$.) The relative complex velocity (2.8) is related to the relative complex potential $F(z)\stackrel {\small {\rm def}}= \varPhi (x,y) + \mathrm {i}\varPsi (x,y)$ by $W=\mathrm {d} F/\mathrm {d} z$, where $U=\varPhi _x=\varPsi _y$ and $V=\varPhi _y=-\varPsi _x$ (see Clamond et al. (Reference Clamond, Labarbe and Henry2023) for more details).

Following Clamond & Constantin (Reference Clamond and Constantin2013), we introduce a complex ‘pressure’ function as

(2.9)\begin{equation} {\mathfrak P}(z) \stackrel{\text{def}}= gd + \tfrac{1}{2}{B}_{s} + \omega ({\psi}_{s} - {\psi}_{b}) - \tfrac{1}{2} W(z)^2 = gd + \tfrac{1}{2}{B}_{b} - \tfrac{1}{2} W(z)^2, \end{equation}

that is holomorphic in the fluid domain $\varOmega$, its restriction to the flat bed $y=-d$ having zero imaginary part and real part $p_{b}$. Thus, $p_{b}$ determines $\mathfrak {P}$ uniquely throughout the fluid domain, i.e. ${\mathfrak P}(z)={p}_{b}(z+\textrm{i} d )$. Note that $p$ introduced in (2.2a,b) coincides with the real part of $\mathfrak {P}$ only on $y=-d$ because the former is not a harmonic function in the fluid domain (Constantin Reference Constantin2006; Constantin & Strauss Reference Constantin and Strauss2010).

As for irrotational waves, it is useful (Clamond Reference Clamond2013, Reference Clamond2018; Clamond et al. Reference Clamond, Labarbe and Henry2023) to introduce the anti-derivative of $\mathfrak {P}(z)$,

(2.10)\begin{equation} \mathfrak{Q}(z) \stackrel{\text{def}}= \int_{z_0}^z \left[ \mathfrak{P}(z') - gd \right] \textrm{d} z' =\int_{z_0}^z \tfrac{1}{2} \big[ {B}_{b} - W(z')^2 \big] \textrm{d} z', \end{equation}

where $z_0$ is an arbitrary constant. For the same abscissa $x$, the functions $\mathfrak {Q}$ at the free surface (i.e. ${\mathfrak {Q}}_{s}(x)$) and at the bottom (i.e. ${\mathfrak {Q}}_{b}(x)$) satisfy the relation

(2.11)\begin{equation} {\mathfrak{Q}}_{s}(x) - {\mathfrak{Q}}_{b}(x) = \int_{x-\textrm{i} d}^{x+\textrm{i}\eta(x)} \left[ \mathfrak{P}(z)- gd \right] \textrm{d} z = \frac{\textrm{i} h(x) {B}_{b}}{2} - \int_{x-\textrm{i} d}^{x+\textrm{i}\eta(x)} \frac{W(z)^2}{2} \textrm{d} z. \end{equation}

2.3. Cauchy integral formula

In the complex $z$ plane the boundaries are analytic curves defined by ${z}_{s}=x+\mathrm {i}\eta$ and ${z}_{b}=x-\mathrm {i} d$. For a holomorphic function $\varXi (z)$, the Cauchy integral formula applied to the fluid domain $\varOmega$ (assuming non-intersecting and non-overturning seabed and free surface) yields

(2.12)\begin{equation} \textrm{i}\vartheta \varXi(z) = \text{P.V.} \oint \frac{\varXi(z')}{z'-z} \textrm{d} z' = \int_{-\infty}^\infty \frac{{\varXi}_{b}'\,\textrm{d}\,x'}{{z}_{b}'-z} - \int_{-\infty}^\infty \frac{(1 + \textrm{i}\eta_x') {\varXi}_{s}'\,\textrm{d}\,x'}{{z}_{s}'-z}, \end{equation}

with primes denoting the dependence on the dummy variable (e.g. ${\varXi }_{b}'\stackrel {\small {\rm def}}={\varXi }_{b}(x')$) and where $\vartheta =\{2{\rm \pi},0,{\rm \pi} \}$ respectively inside, outside and at the smooth boundary of the domain. We emphasis that, in this paper all integrals must be taken in the sense of the Cauchy principal value (P.V.), even if it is not explicitly mentioned for brevity. When $\operatorname {Im}\{{\varXi }_{b}\}=0$, the bottom boundary condition can be taken into account with the method of images, i.e. exploiting the Schwarz reflection principle (Morse & Feshbach Reference Morse and Feshbach1953), yielding

(2.13)\begin{equation} \varXi(z) = \frac{\textrm{i}}{\vartheta} \int_{-\infty}^{\infty} \frac{(1+\textrm{i}\eta_x') {\varXi}_{s}'\,\textrm{d}\,x' }{ {z}_{s}' - z } - \frac{\textrm{i}}{\vartheta} \int_{-\infty}^{\infty} \frac{(1- \textrm{i}\eta_x'){\bar{\varXi}}_{s}'\,\textrm{d}\,x' }{ {\bar{z}}_{s}' - z - 2\textrm{i} d}, \end{equation}

where an overbar denotes the complex conjugation. Note that (2.13) is valid in finite depth (provided that $\operatorname {Im}\{{\varXi }_{b}\}=0$), and in infinite depth if ${\varXi }_{b}\to 0$ as $d\to \infty$. Examples of functions satisfying these conditions are $\varXi =W+c_1$, $\varXi =W^2-{B}_{b}$ and $\varXi = \mathfrak {P}-gd$, so (2.13) provides an expression for computing these functions in arbitrary depth.

2.4. Integral formulations for periodic waves

For $L$-periodic waves (with $L=2{\rm \pi} /k$), the Cauchy kernel is periodised in the horizontal direction, along the interval $\mathcal {I}=[0,L]$, leading to the Cauchy integral formula with Hilbert kernel (Gakhov Reference Gakhov1990)

(2.14)\begin{equation} \varXi(z) = \frac{\textrm{i} k}{2\vartheta} \int_{\mathcal{I}} \left[\cot\!\left( k\frac{{z}_{s}'-z}{2} \right) (1+\textrm{i}\eta_x') {\varXi}_{s}' - \cot\!\left( k\frac{{z}_{b}'-z}{2} \right){\varXi}_{b}' \right] \textrm{d}\,x' . \end{equation}

Alternatively, using the method of images, along with the identity (A10), the Cauchy integral can be written as

(2.15)\begin{align} \varXi(z) &= \frac{k}{\vartheta} \int_{\mathcal{I}} \left[ \mbox{Li}_{0} \{\exp({\textrm{i} k ({z}_{s}'-z)})\} (1+\textrm{i}\eta_x') {\varXi}_{s}' \right. \nonumber\\ &\quad \left.+ \mbox{Li}_{0}\{\exp({\textrm{i} k (z-{\bar{z}}_{s}'+2\textrm{i} d)})\} (1-\textrm{i}\eta_x') {\bar{\varXi}}_{s} \right] \textrm{d}x' \nonumber\\ &\quad + {\rm \pi}\vartheta^{{-}1} \langle (1+\textrm{i}\eta_x){\varXi}_{s} + (1-\textrm{i}\eta_x){\bar{\varXi}}_{s} \rangle, \end{align}

where $\mbox{Li}_{\nu }$ is the $\nu$th polylogarithm whose definition is given in Appendix A along with useful relations. It is worth noting that the last term on the right-hand side of (2.15) corresponds to the zeroth Fourier coefficient (not present when the holomorphic function $\varXi (z)$ has zero mean over the wave period). Equation (2.15) can be rewritten using the identity (A9), yielding

(2.16)\begin{align} \varXi(z) &= \frac{\textrm{i}}{\vartheta} \int_{\mathcal{I}} \frac{\partial}{\partial z} \left[\mbox{Li}_{1} \{\exp({\textrm{i} k ({z}_{s}'-z)})\} (1+\textrm{i}\eta_x') {\varXi}_{s}' \right. \nonumber\\ &\quad \left. -\mbox{Li}_{1} \{\exp({\textrm{i} k (z-{ \bar{z}}_{s}'+2\textrm{i} d)})\} (1-\textrm{i}\eta_x') {\bar{\varXi}}_{s}\right] \textrm{d}\,x' \nonumber\\ &\quad + 2{\rm \pi}\vartheta^{{-}1} \left\langle \operatorname{Re}{\left\{{\varXi}_{s}\right\}} - \eta_x \operatorname{Im}{\left\{{\varXi}_{s}\right\}} \right\rangle. \end{align}

At the free surface (where $\vartheta ={\rm \pi}$, $z={z}_{s}$ and $\mathrm {d} {z}_{s} = (1+\mathrm {i}\eta _x)\,\mathrm {d}\kern 0.06em x$), carefully applying the Leibniz integral rule (c.f. (B4) in Appendix B) on the singular term in the integrand, (2.16) reduces to

(2.17)\begin{align} \left(1+\textrm{i}\eta_x\right){\varXi}_{s} &= \frac{\textrm{i}}{2{\rm \pi}} \frac{\textrm{d}}{\textrm{d}\,x} \int_{\mathcal{I}} \left[ \mbox{Li}_{1} \{ \exp({\textrm{i} k ({z}_{s}'-{z}_{s})}) \} (1+\textrm{i}\eta_x') {\varXi}_{s}' \right. \nonumber\\ &\quad \left. -\mbox{Li}_{1} \{ \exp({\textrm{i} k ({z}_{s}-{\bar{z}}_{s}'+2\textrm{i} d)}) \} (1-\textrm{i}\eta_x') {\bar{\varXi}}_{s} \right] \textrm{d}\,x' \nonumber\\ &\quad + \left(1+\textrm{i}\eta_x\right) \left\langle \operatorname{Re}{\left\{{\varXi}_{s}\right\}} - \eta_x \operatorname{Im}{\left\{{\varXi}_{s}\right\}} \right\rangle . \end{align}

It should be noted that, obviously, aperiodic equations can be obtained from the periodic ones letting $L\to \infty$ (i.e. $k\to 0^+$). Thus, from now on, we only consider periodic waves.

4.1. Eulerian formulation

From expressions (2.4) and (2.8), the complex (irrotational part of the) velocity at the surface is given explicitly by

(4.1)\begin{equation} {W}_{s} = \omega h + \sigma \left(1-\textrm{i}\eta_x\right)\sqrt{({B}_{s}-2g\eta)/(1+\eta_x^2)}, \end{equation}

with $\sigma \stackrel {\small {\rm def}}=\mp 1$ denoting waves propagating upstream or downstream, respectively. The parameter $\sigma$ is introduced for convenience in order to characterise the (arbitrarily chosen) direction of the wave propagation in a ‘fixed’ frame of reference, i.e. $\sigma =-1$ if the wave travels toward the increasing $x$ direction in this frame and, obviously, $\sigma =+1$ if the wave travels toward the decreasing $x$ direction.

Considering the holomorphic function $\varXi =W+c$ ($c$ being an arbitrary definition of the phase speed), the left-hand side of (2.17) follows directly from (4.1),

(4.2)\begin{equation} \left(1+\textrm{i}\eta_x\right){\varXi}_{s} = \left(\omega h + c\right) \left(1+\textrm{i}\eta_x\right) + \sigma \sqrt{({B}_{s}-2g\eta)(1+\eta_x^2)}, \end{equation}

where the radicand is purely real since ${B}_{s}\geqslant \max (2g\eta )$ for all waves. Substituting expression (4.2) in (2.17), the integral term splits into several contributions as

(4.3) \begin{align} \omega \eta \left(1+\textrm{i}\eta_x\right) &= \frac{\textrm{i}}{2{\rm \pi}}\frac{\textrm{d}}{\textrm{d}\,x} \left[\int_{\mathcal{C}} \mbox{Li}_{1}(\exp({\textrm{i} k ({z}_{s}'-{z}_{s})})) (\omega h' + c) \textrm{d}{z}_{s}' \right. \nonumber\\ &\quad + \left. \int_{\mathcal{C}}\mbox{Li}_{1}(\exp({\textrm{i} k ({z}_{s}-{\bar{z}}_{s}'+2\textrm{i} d)})) (\omega h' + c) \textrm{d} {\bar{z}}_{s}' \right. \nonumber\\ &\quad + \left. \sigma\int_{\mathcal{I}} \{ \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}'-{z}_{s})})) - \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}-{\bar{z}}_{s}'+2\textrm{i} d)}))\} \right. \nonumber\\ &\quad \vphantom{\int_{\mathcal{C}}} \times \left. \sqrt{({B}_{s}-2g\eta')(1+\eta_x'^2)}\,\textrm{d}\,x' \right] \nonumber\\ &\quad + \left(1+\textrm{i}\eta_x\right) \langle {u}_{s} (1+\eta_x^2) \rangle - \sigma \sqrt{({B}_{s}-2g\eta)(1+\eta_x^2)}, \end{align}

where $\mathcal {C}$ represents the free surface path, i.e. we use the brief notations $\int _{\mathcal {C}}(\cdots )\,\mathrm {d}{z}_{s}'\stackrel {\small {\rm def}}= \int _{\mathcal {I}}(\cdots )(1+\mathrm {i}\eta _x')\,\mathrm {d}\kern 0.06em x'$ and $\int _{\mathcal {I}}(\cdots )\,\mathrm {d}\kern 0.06em x'\stackrel {\small {\rm def}}=\int _0^L(\cdots )\,\mathrm {d}\kern 0.06em x'$. The first term inside the square bracket on the right-hand side of (4.3) reduces to

(4.4)\begin{align} J_1 &\stackrel{\text{def}}= \int_{\mathcal{C}} \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}'- {z}_{s})})) (\omega h' + c)\,\textrm{d}{z}_{s}' = \omega \int_{\mathcal{C}} \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}'- {z}_{s})})) \eta'\,\textrm{d}{z}_{s}' \nonumber\\ &= \omega \int_{\mathcal{I}} \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}'- {z}_{s})}) ) \eta' (1+\textrm{i}\eta_x') \textrm{d}\,x'= \frac{\textrm{i}\omega}{k} \int_{\mathcal{I}} \mbox{Li}_{2}(\exp({\textrm{i} k({z}_{s}'-{z}_{s})}))\eta_x'\,\textrm{d}\,x' \nonumber\\ &={-} \frac{\omega}{k} \int_{\mathcal{I}} \mbox{Li}_{2}(\exp({\textrm{i} k({z}_{s}'-{z}_{s})}) ) \textrm{d}\,x', \end{align}

where we exploited the property (Clamond Reference Clamond2018)

(4.5)\begin{equation} \int_{\mathcal{C}} \mbox{Li}_{\nu}\!\left( \exp({\textrm{i} k{z}_{s}}) \right)\textrm{d} {z}_{s} = 0. \end{equation}

Similarly, we have

(4.6)\begin{align} J_2 &\stackrel{\text{def}}= \int_{\mathcal{C}} \mbox{Li}_{1}(\exp({\textrm{i} k({z}_{s}- {\bar{z}}_{s}'+2\textrm{i} d)})) (\omega h'+c)\,\textrm{d}{\bar{z}}_{s}' \nonumber\\ &={-}\frac{\omega}{k} \int_{\mathcal{I}} \mbox{Li}_{2}(\exp({\textrm{i} k({z}_{s}-{\bar{z}}_{s}'+2\textrm{i}d)})) \textrm{d}\,x'. \end{align}

Finally, after integrating the whole expression (4.3) and retaining the imaginary part only, we obtain the equation for the free surface

(4.7)\begin{align} K &= \frac{\omega\eta^2}{2} - \eta \langle {u}_{s} (1+\eta_x^2) \rangle + \frac{\omega}{2{\rm \pi} k} \int_{\mathcal{I}} \operatorname{Re}\{\mathcal{L}_{2}\}\, \textrm{d}\,x' \nonumber\\ &\quad - \frac{\sigma}{2{\rm \pi}} \int_{\mathcal{I}} \operatorname{Re}\{\mathcal{L}_{1}\} \sqrt{( {B}_{s}-2g\eta') (1+\eta_x'^2)}\, \textrm{d}\,x', \end{align}

where $K$ is an integration constant obtained enforcing the mean-level condition (2.1), i.e.

(4.8)\begin{equation} K = \left\langle\frac{\omega\eta^2}{2} + \frac{\omega}{2{\rm \pi} k} \int_{\mathcal{I}} \operatorname{Re}\{\mathcal{L}_{2}\}\,\textrm{d}\,x'- \frac{\sigma}{2{\rm \pi}} \int_{\mathcal{I}} \operatorname{Re}\{\mathcal{L}_{1}\} \sqrt{( {B}_{s}-2g\eta') (1+\eta_x'^2)}\,\textrm{d}\,x'\right\rangle. \end{equation}

From now on, the same notation $K$ is used to denote integration constants in the different surface recovery formulae. The exact values of these constants is obtained, as above, by enforcing the condition (2.1). Moreover, we have introduced, for brevity, the notation for the kernels

(4.9)\begin{equation} \mathcal{L}_{\nu} \stackrel{\text{def}}= \mbox{Li}_{\nu}(\exp({\textrm{i} k({z}_{s}'- {z}_{s})})) - \mbox{Li}_{\nu}(\exp({\textrm{i} k({z}_{s}- {\bar{z}}_{s}'+2\textrm{i}d)})). \end{equation}

Equation (4.7) is a nonlinear integro-differential equation for the computation of the surface elevation $\eta$. Once $\eta$ is obtained solving (4.7) numerically, one can compute the corresponding bottom pressure. This bottom pressure can then be considered as ‘experimental data’ to illustrate the reconstruction procedure described below. This is necessary when such data are not available from physical measurements, as is the case with rotational waves of arbitrary shape. We can only stress how valuable in situ measurements would be to test the recovery of surface wave profiles from realistic data.

4.2. Lagrangian formulation

Expression (3.11) provides an implicit definition of $T_L$ if the wavelength $L$ is fixed a priori. It is now of interest to convert our formula for the surface elevation (4.7) using the Lagrangian description introduced previously (the variable of interest now becomes ${\theta }_{s}$). After some simplifications, we obtain

(4.10)\begin{align} & \frac{\omega\sigma}{2} \eta^2 - \frac{k\eta}{2{\rm \pi}} \int_{0}^{T_L} ({B}_{s} - 2g\eta)\,\textrm{d} t - \frac{1}{2{\rm \pi}} \int_{0}^{T_L} \operatorname{Re}{\{\mathcal{L}_1\}} ({B}_{s} - 2g\eta')\,\textrm{d} t' \nonumber\\ &\quad + \frac{\omega\sigma}{2{\rm \pi} k} \int_{0}^{T_L} \operatorname{Re}{\{\mathcal{L}_2\}} \cos({{\theta}_{s}'}) \sqrt{{B}_{s} - 2g\eta'}\, \textrm{d} t' = K, \end{align}

where $K$ is recovered using condition (2.1) in the same way that (4.8) was obtained.

The computation of (4.10) involves weak logarithmic singularities in the kernel of the $\mathcal {L}_1$ operator. In the numerical implementation we use a similar approach as the one presented by Clamond (Reference Clamond2018), by subtracting the regular part of the operator. Thence, we obtain an explicit expression for the regularised finite integral, i.e.

(4.11)\begin{align} \int_{0}^{T_L} \operatorname{Re}{\{\mathcal{L}_1\}} ({B}_{s} - 2g\eta')\,\textrm{d} t' &={-}2g \int_{0}^{T_L} \operatorname{Re}{\{\mathcal{L}_1\}} (\eta' - \eta)\,\textrm{d} t' \nonumber\\ &\quad +({B}_{s} - 2g\eta) \int_{0}^{T_L} \operatorname{Re}{\{\mathcal{L}_1 - \hat{\mathcal{L}}_1\}}\,\textrm{d} t', \end{align}

where we introduced

(4.12)\begin{equation} \hat{\mathcal{L}}_{\nu} \stackrel{\text{def}}= \mbox{Li}_{\nu} \{\exp({\textrm{i} \tau(t'-t)})\} - \mbox{Li}_{\nu}\{ \exp({\textrm{i} \tau(t-t')}) \exp({-2kd}) \}, \end{equation}

with $\tau \stackrel{\text{def}}{=} 2{\rm \pi} /T_L$. Considering $t\to t'$ in both integrands in (4.11), we have

(4.13)$$\begin{gather} \lim_{t\to t'} \operatorname{Re}{\{\mathcal{L}_1\}} (\eta - \eta') = 0 , \end{gather}$$

(4.14)$$\begin{gather}\lim_{t\to t'} \operatorname{Re}{\{\mathcal{L}_1 - \hat{\mathcal{L}}_1\}} = \log\!\left[\left(\frac{1 - \exp({-2kh})}{1 - \exp({-2kd})} \right) \frac{\tau/k}{\sqrt{{B}_{s} - 2g\eta}}\right]. \end{gather}$$

The numerical computation of Lagrangian surface waves is thus done by solving the nonlinear expression (4.10), providing a given wave amplitude $\eta |_{t=0} = a$. The wave period $T_L$, as well as the Bernoulli constant $B_s$, are both obtained from the Lagrangian counterpart of (2.1) and the requirement that $\int _{0}^{T_L} \mathrm {d} {z}_{s} = 2{\rm \pi} \sigma /k$. Moreover, since the abscissa $x(t)$ is given explicitly from the definition of the wave profile $\eta (t)$ through the real part of (3.7), we only have to solve expression (4.10) for $\eta$ – and not for both $x$ and $\eta$, as done previously by Da Silva & Peregrine (Reference Da Silva and Peregrine1988) and by Vanden-Broeck (Reference Vanden-Broeck1994) – thus allowing us to further reduce the numerical cost.

5.1. Eulerian boundary integral formulation

Let us consider the Cauchy integral formula (2.17) (without the method of images) for the holomorphic function $\varXi (z) = \mathfrak {P}(z) - gd$. It yields

(5.1)\begin{align} \mathfrak{P} - gd &= \frac{k}{\vartheta} \int_{\mathcal{I}} \left[\mbox{Li}_{0}(\exp({\textrm{i} k ({z}_{s}'-z)})) (1+\textrm{i}\eta_x') ({\mathfrak{P}}_{s}' - gd) \right. \nonumber\\ &\quad \left.- \mbox{Li}_{0}(\exp({\textrm{i} k ({z}_{b}'-z)})) ({\mathfrak{P}}_{b}'- gd) \right] \textrm{d}\,x' \nonumber\\ &\quad + \langle (1+\textrm{i}\eta_x) ({\mathfrak{P}}_{s} - gd) - ({\mathfrak{P}}_{b}' - gd) \rangle . \end{align}

In order to simplify the latter expression, we exploit both definitions of ${\mathfrak {P}}_{b}\stackrel {\small {\rm def}}= p_b(x)$ and $\mathrm {d} {\mathfrak {Q}}_{s}/\mathrm {d}\kern 0.06em x = ({\mathfrak {P}}_{s}-gd)(1+ \mathrm {i}\eta _x)$ (we recall that ${\mathfrak {Q}}_{s}$ is a periodic function). Hence, expression (5.1) can be rewritten as

(5.2)\begin{align} \mathfrak{P} - gd &= \frac{\textrm{i}}{\vartheta} \int_{\mathcal{I}} \frac{\partial}{\partial z} \left[\mbox{Li}_{1}(\exp({\textrm{i} k ({z}_{s}'-z)})) (1+\textrm{i}\eta_x') ({\mathfrak{P}}_{s}' - gd)\right. \nonumber\\ &\quad \left. -\mbox{Li}_{1}(\exp({\textrm{i} k ({z}_{b}'-z)})) ({p}_{b}' - gd)\right] \textrm{d}\,x'. \end{align}

Before continuing, let the compact notations for the polylogarithmic kernels be

(5.3a,b)\begin{equation} \mathcal{K}_{\nu} \stackrel{\text{def}}= \mbox{Li}_{\nu}( \exp({\textrm{i} k ({z}_{s}'-{z}_{s})})) \quad \textrm{and} \quad \mathcal{J}_{\nu} \stackrel{\text{def}}= \mbox{Li}_{\nu}( \exp({\textrm{i} k ({z}_{b}'-{z}_{s})})) . \end{equation}

We now evaluate expression (5.2) at the free surface, reducing it to

(5.4)\begin{equation} ({\mathfrak{P}}_{s} - gd) (1+\textrm{i}\eta_x) = \frac{\textrm{i}}{2{\rm \pi}} \frac{\textrm{d}}{\textrm{d}\,x} \int_{\mathcal{I}} \left[ \mathcal{K}_{1} (1+\textrm{i}\eta_x') ({\mathfrak{P}}_{s}' - gd) - \mathcal{J}_{1} ({\mathfrak{P}}_{b}' - gd) \right] \textrm{d}\,x'. \end{equation}

As one can notice, the left-hand side of (5.4) is the integrand of ${\mathfrak {Q}}_{s}$. For that reason, we first integrate the whole expression over the $x$ coordinate and then split the contributions from the surface and the bottom,

(5.5)\begin{equation} {\mathfrak{Q}}_{s} = \frac{\textrm{i}}{2{\rm \pi}} \int_{\mathcal{I}} \mathcal{K}_{1} (1+\textrm{i}\eta_x') ({\mathfrak{P}}_{s}' - gd) \,\textrm{d}\,x' - \frac{\textrm{i}}{2{\rm \pi}} \int_{\mathcal{I}} \mathcal{J}_{1} ({p}_{b}' - gd) \,\textrm{d}\,x' + K, \end{equation}

for a constant of integration $K$ also obtained by applying condition (2.1).

The next step is to decompose the integrand in the first integral of (5.5). To do so, we merely replace the complex velocity by its expression (2.8) and exploit the identity (4.5) to cancel out some constant terms, i.e.

(5.6)\begin{align} & \int_{\mathcal{C}} \mathcal{K}_{1} ({\mathfrak{P}}_{s}' - gd) \,\textrm{d} {z}_{s}' \nonumber\\ &\quad ={-} \frac{1}{2}\int_{\mathcal{C}} \mathcal{K}_{1} {W}_{s}'^2 \,\textrm{d} {z}_{s}' ={-} \frac{1}{2} \int_{\mathcal{C}}\mathcal{K}_{1}\left[ \omega h' + (1-\textrm{i}\eta_x') \sqrt{\frac{({B}_{s}-2g\eta')}{(1+\eta_x'^2)}} \right]^2 \textrm{d} {z}_{s}' \nonumber\\ &\quad ={-}\frac{\omega^2}{2} \int_{\mathcal{C}}\mathcal{K}_{1} h'^2 \,\textrm{d} {z}_{s}' - \omega\sigma \int_{\mathcal{C}} \mathcal{K}_{1} h' (1-\textrm{i}\eta_x') \sqrt{\frac{({B}_{s}-2g\eta')}{(1+\eta_x'^2)}}\,\textrm{d} {z}_{s}' \nonumber\\ &\qquad -\frac{1}{2} \int_{\mathcal{C}} \mathcal{K}_{1} (1-\textrm{i}\eta_x')^2 \frac{({B}_{s}-2g\eta')}{(1+\eta_x'^2)} \textrm{d} {z}_{s}'. \end{align}

After some algebraic manipulations, it further reduces to

(5.7)\begin{align} \int_{\mathcal{C}} \mathcal{K}_{1} ({\mathfrak{P}}_{s}' - gd)\,\textrm{d} {z}_{s}' &={-}\frac{\omega^2}{2} \int_{\mathcal{C}} \mathcal{K}_{1} \eta'^2 \,\textrm{d}{z}_{s}'- \omega^2d \int_{\mathcal{C}} \mathcal{K}_{1} \eta' \,\textrm{d}{z}_{s}' \nonumber\\ &\quad -\frac{1}{2} \int_{\mathcal{I}} \mathcal{K}_{1} ({B}_{s} - 2g\eta') (1 - \textrm{i}\eta_x') \,\textrm{d}\,x' \nonumber\\ &\quad -\omega\sigma \int_{\mathcal{I}} \mathcal{K}_{1} h' \sqrt{({B}_{s}- 2g\eta')(1+\eta_x'^2)} \,\textrm{d}\,x' \nonumber\\ &={-}\frac{\textrm{i}\omega^2}{k} \int_{\mathcal{I}} \mathcal{K}_{2} \eta' \eta_x' \,\textrm{d}\,x' + \frac{\omega^2d + g}{k} \int_{\mathcal{I}} \mathcal{K}_{2} \,\textrm{d}\,x' \nonumber\\ &\quad - \int_{\mathcal{I}} \mathcal{K}_{1} ({B}_{s} - 2g\eta') \,\textrm{d}\,x' \nonumber\\ &\quad -\omega\sigma \int_{\mathcal{I}} \mathcal{K}_{1} h' \sqrt{({B}_{s}- 2g\eta')(1+\eta_x'^2)} \,\textrm{d}\,x'. \end{align}

Substituting (5.7) into (5.5), the imaginary part yields the Eulerian integral formulation for the surface recovery

(5.8)\begin{align} 2{\rm \pi} \operatorname{Im} \{ {\mathfrak{Q}}_{s} \} &= \frac{\omega^2}{2k} \int_{\mathcal{I}} \operatorname{Im} \{ \mathcal{K}_{2} \} (\eta'^2)_x \,\textrm{d}\,x' - \int_{\mathcal{I}} \operatorname{Re} \{ \mathcal{K}_{1} \} ({B}_{s} - 2g\eta') \,\textrm{d}\,x' \nonumber\\ &\quad +\frac{\omega^2d + g}{k} \int_{\mathcal{I}} \operatorname{Re} \{ \mathcal{K}_{2} \}\,\textrm{d}\,x'- \int_{\mathcal{I}} \operatorname{Re} \{ \mathcal{J}_{1} \} ({\mathfrak{P}}_{b}' - gd) \,\textrm{d}\,x' \nonumber\\ &\quad - \omega\sigma \int_{\mathcal{I}} \operatorname{Re} \{ \mathcal{K}_{1} \} h' \sqrt{({B}_{s} - 2g\eta')(1+\eta_x'^2)}\,\textrm{d}\,x' - 2{\rm \pi}\operatorname{Im}\{K\}, \end{align}

where $K$ is the same constant as in (5.5).

Equation (5.8) is a nonlinear integral equation for the free surface recovery from the bottom pressure. Being strictly Eulerian, this equation is not suitable for overhanging waves. For the latter, one can proceed as detailed in the following section.

5.2. Hybrid formulation

Equation (5.8) involves integrals at the free surface and at the bottom. In practice, the bottom pressure is given at some known abscissa $x$, so the bottom integral must be kept in Eulerian form. However, Eulerian integrals are not suitable for overhanging waves, so we rewrite surface integrals in their Lagrangian counterparts. Doing so, the inverse problem is described by a mixed Eulerian–Lagrangian formalism.

Thus, by rewriting the surface integrals in (5.8) in the Lagrangian description and keeping the bottom integral in Eulerian form, one gets the general expression for the surface recovery:

(5.9)\begin{align} 2{\rm \pi} \operatorname{Im} \{ {\mathfrak{Q}}_{s} \} &= \frac{\omega^2}{k} \int_{0}^{T_L} \left[k^{{-}1}\operatorname{Re} \{ \mathcal{K}_{3} \} + \left( h' + g\omega^{{-}2} \right) \operatorname{Re} \{ \mathcal{K}_{2} \}\right] \cos({\theta}_{s}') \sqrt{{B}_{s} - 2g\eta'} \,\textrm{d} t' \nonumber\\ &\quad - \int_{0}^{T_L} \operatorname{Re} \{ \mathcal{K}_{1} \} \left[ \cos({\theta}_{s}') \sqrt{{B}_{s} - 2g\eta'} + \sigma\omega h' \right]\left({B}_{s} - 2g\eta'\right)\textrm{d} t' \nonumber\\ &\quad - \int_{0}^L \operatorname{Re} \{ \mathcal{J}_{1} \} ({p}_{b}' - gd)\,\textrm{d}\,x' - 2{\rm \pi}\operatorname{Im}\{K\} . \end{align}

This reformulation is necessary for practical recovery of overhanging waves. It is, of course, also suitable for non-overhanging waves.

6.1. Details on the overall methodology

From a numerical standpoint, the nonlinear integral equations (4.7) and (5.9), respectively employed for computing the surface profile and recovering it from the bottom pressure, possess notable characteristics. First, these equations eliminate the need for evaluating the derivative of ${\theta }_{s}$ at any point. Second, we can rely on Fourier analysis since the kernels are periodic and use the trapezoidal rule for numerical integration.

Since we already know the value of $g$ and the wavenumber $k$ (or equivalently, the wavelength $L$) as inputs in our numerical scheme, we can deduce the depth of the layer $d$ from the definition of the hydrostatic law (Clamond Reference Clamond2013; Clamond et al. Reference Clamond, Labarbe and Henry2023). For simplicity, we consider that the constant vorticity $\omega$ is known. However, $\omega$ can also be determined by adapting the procedure described by Clamond et al. (Reference Clamond, Labarbe and Henry2023). Then, initialising our numerical schemes with an appropriate initial guess, either from linear theory or by a previous iteration, as done by Da Silva & Peregrine (Reference Da Silva and Peregrine1988), we observe fast convergence for a discrete set of $N$ equidistant points. In our simulations we typically use $N=128$ for moderate waves and $N=512$ for waves (regardless of whether they exhibit overturning or not) with large amplitudes.

When investigating overhanging waves, the non-algebraic nonlinearities inherent in the problem often pose numerical challenges. The most troublesome issue arises from aliasing errors present in the functions spectra, a phenomenon occasionally referred to as ‘spectral blocking’ (Boyd Reference Boyd2001), leading to exponential growth of high frequencies. As a consequence, this aliasing effect prevents the spectral accuracy inherent in our formulation. To address this issue, we employ the ‘zero padding’ method, which involves increasing the size of the quadrature in Fourier space while appending zeros above the Nyquist frequency. Subsequently, we transform the functions back to physical space, compute the nonlinear terms and filter out the previously introduced zero frequencies from the spectrum. For quadratic nonlinear terms, enlarging the degree of quadrature by a factor of $3/2$ has proven sufficient to mitigate this phenomenon (Patterson & Orszag Reference Patterson and Orszag1971). However, in the context of non-algebraic nonlinearities encountered in this study, this argument does not hold, and the exact value of the enlargement factor remains unknown. Instead, we utilize a factor of $2$ (typically suitable for cubic nonlinearities) when performing products throughout the algorithm. Although we experimented with larger factors in our simulations, it appeared that this value was adequate for eliminating spurious frequencies in most aliased spectra.

In addition to aliasing, we noticed that it is numerically more efficient to perform a change of coordinates in evaluating the finite integrals. Particularly, significant variations in ${\theta }_{s}$ occur where the wave undergoes overturning, indicating a requirement for additional collocation points in these regions. Thus, rather than evaluating the previous integrals with respect to the time variable $t\in [0,T_L]$, we introduce a new integration variable $\varXi (t)$, defined as

(6.1)\begin{equation} \frac{\textrm{d} \xi}{\textrm{d} t} \stackrel{\text{def}}= \sqrt{1 + \beta\sin\!\left({\theta}_{s}\right)^2}, \end{equation}

where $\beta$ is a scaling parameter arbitrarily set. A similar change of coordinate can be found in Da Silva & Peregrine (Reference Da Silva and Peregrine1988). In most simulations involving overhanging waves, we employ a value of $\beta =2{\rm \pi} /(kdN)$.

Finally, we solve the whole system of (4.7) and (5.9) with the built-in iterative solver fsolve from the Matlab software, using the Levenberg–Marquardt algorithm.

6.2. Rotational periodic overhanging wave

The first case of interest, as shown in figure 2, depicts a periodic overturning wave of large amplitude. This particular scenario poses significant computational challenges, making it an ideal benchmark for evaluating the robustness of our recovery procedure. Utilizing the pressure data highlighted in figure 2(a), we successfully reconstruct the surface profile using the expression 5.9, resulting in excellent agreement, as evident from figure 2(b). Notably, the change of variables achieved through (6.1) effectively concentrates the points in regions where $\theta _s$ varies the most. Consequently, we attain a high level of accuracy, with $\vert \vert \eta - \eta ^{ex}\vert \vert _{\infty } \approx 4.48\times 10^{-3}$, where $\eta ^{ex}$ corresponds to a numerical solution obtained by solving (4.10). Furthermore, for the computation of the unknown Bernoulli constant at the surface, the error is approximately $\vert {B}_{s} - {B}_{s}^{ex} \vert \approx 5.83\times 10^{-3}$.

Figure 2. Numerical demonstration of the surface recovery procedure for an overhanging and periodic steady wave with a constant vorticity $\omega \sqrt {d/g}=3\sqrt {2}$. (a) Eulerian representation of the bottom pressure. The dashed-dot line corresponds to the mean bottom pressure. (b) Surface wave profile (blue circles) obtained from expression (4.7). The red line represents the surface reconstruction achieved from ${p}_{b}$ through (5.9), while the dashed-dot line indicates the mean water level.

We emphasise the effectiveness of our recovery procedure by presenting the velocity field within the fluid layer in both upper panels of figure 3. Notably, a pair of stagnation points is observed in figure 3(c) on the flat bottom boundary, which often poses challenges when utilizing conformal mapping techniques. In the present context, since our methodology operates solely in the physical plane, our general approach can readily reconstruct the surface profile regardless of the presence or absence of stagnation points.

Figure 3. Isovalues of (a) $u/\sqrt {gd}$, (b) $v/\sqrt {gd}$, (c) $\log {(\vert w\vert /\sqrt {gd})}$ and (d) $\arg {( w/\sqrt {gd})}$ for the configuration displayed in figure 2.

6.3. Rotational solitary surface wave

Computing solitary waves using our procedure is straightforward but necessitates considering a relatively large numerical domain to accurately capture their behaviour far from the crest. Meanwhile, to ensure that the wave is located above the mean water level, we substitute the condition (2.1) with

(6.2)\begin{equation} \eta^{(n+1)} = \eta^{(n)} - \min\eta^{(n)}, \end{equation}

instead, at each numerical iteration $n$.

The recovery process for a solitary wave is presented in figure 4, which illustrates the pressure data in the upper panel and the corresponding surface profile in the lower panel. We consider here a solitary wave of relatively small amplitude because it is quite challenging. Indeed, the larger the solitary wave, the faster its decay (thus requiring a smaller computational box) and the larger the ratio signal/noise (for field data). So, in this respect, the recovery of small amplitude solitary waves is more challenging.

Figure 4. Surface recovery procedure for a solitary steady wave with a constant vorticity $\omega \sqrt {d/g}=1$. Both panels have the same legend as figure 2.

In the present specific case, the agreement with a given numerical solution is great, yielding numerical errors of $\vert \vert \eta - \eta ^{ex} \vert \vert _{\infty } \approx 1.46\times 10^{-4}$ and $\vert {B}_{s} - {B}_{s}^{ex}\vert \approx 4.21\times 10^{-5}$ for the surface profile and the Bernoulli constant, respectively. Unfortunately, spurious infinitesimal oscillations (located far from the wave crest) prevent us from reaching the much better accuracy that we would expect from our method. In order to facilitate these computations and remove unwanted oscillations, another change of variables can be implemented (similar to the approach (6.1) used for the periodic case) to concentrate the quadrature points near the crest, rather than far away where the elevation is infinitesimal. However, we reserve this task for future investigations, which will provide more comprehensive details on the efficient computation of solitary waves within this context.