General procedure for free-surface recovery from bottom pressure measurements: application to rotational overhanging waves

Abstract

A novel boundary integral approach for the recovery of overhanging (or not) rotational water waves (with constant vorticity) from pressure measurements at the bottom is presented. The method is based on the Cauchy integral formula and on an Eulerian–Lagrangian formalism to accommodate overturning free surfaces. This approach eliminates the need to introduce a priori a special basis of functions, thus providing a general means of fitting the pressure data and, consequently, recovering the free surface. The effectiveness and accuracy of the method are demonstrated through numerical examples.

1. Introduction

Nonlinear water waves have been extensively studied since the mid-eighteenth century, when Euler introduced his eponymous equation. Since then, surface gravity waves have attracted much attention in their modelling, although scientists quickly realised their inherent complexity. This is one reason why physicists and mathematicians are still interested by the richness of this problem, making it an endless source for research in fluid dynamics. For instance, one crucial concern in environmental and coastal engineering is to accurately measure the surface of the sea for warnings about the formation of large waves near coasts and oceanic routes. One solution to this problem is reconstructing the surface using a discrete set of measurements obtained from submerged pressure transducers (Tsai & Tsai 2009). This approach avoids the limitations of offshore buoy systems, which are susceptible to climatic disasters, located on moving boundaries, and lacking accuracy in wave height estimates (Lin & Yang 2020). Consequently, solving the nonlinear inverse problem associated with water waves is a timely request for building practical engineering apparatus that rely on in situ data.

While the hydrostatic theory was originally used to tackle this problem when first formulated (Bergan, Tørum & Traetteberg 1969; Lee & Wang 1985), it was only recently that nonlinear waves began to be addressed (Oliveras et al. 2012). Some works, such as Constantin (2012), considered conformal mapping to successfully obtain reconstruction formulae. However, the numerical cost associated with solving these implicit relations renders conformal mapping inefficient when dealing with real physical data, i.e. when pressure data are given at known abscissas of the physical plane, not of the conformally mapped one. Actually, introducing some suitable holomorphic functions, it is possible to efficiently solve this nonlinear problem while staying within the physical plane for the recovery procedure (Clamond 2013; Clamond & Constantin 2013). These studies demonstrated the convergence of the reconstruction process and the ability to recover waves of maximum amplitude (Clamond & Henry 2020). Furthermore, the method was adapted to handle cases involving linear shear currents (Clamond, Labarbe & Henry 2023), remarkably recovering the unknown magnitude of vorticity alongside the wave profile and associated parameters. This case is notorious for being challenging, as it allows for the presence of critical layers and/or stagnation points in the fluid domain (Wahlen 2009). Note that the existence of nonlinear water waves with constant vorticity and overturning wave profiles has been demonstrated (Constantin & Varvaruca 2011; Constantin, Strauss & Varvaruca 2016). Note also that, for infinitesimal waves, exact recovery formulae were obtained by Henry & Thomas (2018) for steady waves with general vorticity.

The recovery method studied in Clamond (2013), Clamond & Constantin (2013), Clamond & Henry (2020) and Clamond et al. (2023) has two shortcomings, however. First, it was developed for non-overhanging waves, so it cannot directly address these waves that can occur, in particular, in the presence of constant vorticity. Second, part of this reconstruction procedure is based on analytic continuation of some well-chosen eigenfunctions. If for some waves a ‘good’ choice of eigenfunctions is clear a priori, this is not necessarily the case for complicated wave profiles. By ‘good’ choice, we mean a set of eigenfunctions that provides an accurate representation of the free surface with a minimal amount of modes and, at the same time, that can be easily computed. Even though Fourier series (or integrals) can be used in principle, a large number of eigenfunctions may be required to accurately represent the free surface. Since a surface recovery from bottom pressure is intrinsically ill-conditioned, using a large number of Fourier modes may lead to numerical issues. Another basis should then be preferably employed. For instance, for irrotational long waves propagating in shallow water (i.e. cnoidal waves), the use of Jacobian elliptic functions is effective (Clamond 2013; Clamond & Constantin 2013). However, such an alternative basis is not always easily guessed. Thus, it is desirable to derive a reconstruction procedure independent of a peculiar basis of eigenfunctions. Here, we propose a recovery methodology addressing these two shortcoming.

In this paper we derive a general formulation to address the surface recovery problem using a boundary integral formulation. While a similar approach was described by Da Silva & Peregrine (1988) for computing waves with constant vorticity, to the best of our knowledge it has never been applied in the context of a recovery procedure. The Cauchy integral formula, although singular by definition, proves advantageous from a numerical perspective. Dealing with singular kernels, the integral formulation easily allows for the consideration of arbitrary steady surface waves (periodic, solitary, aperiodic) travelling in a linear shear current, without the need to select a peculiar basis of functions to fit the pressure data. Additionally, this method facilitates the parametrisation of the surface profile, enabling the recovery of overhanging waves with arbitrarily large amplitudes.

Overturning profiles are known to be hard to compute accurately (Vanden-Broeck 1994), which presents a significant challenge due to the ill-posed nature of our problem. Nevertheless, by considering a mixed Eulerian–Lagrangian description at the boundaries, we demonstrate the feasibility of the recovery process. To illustrate the robustness of our method, we present two examples of constant vorticity steady waves: a periodic wave with an overturning surface and a solitary wave. In both scenarios, we achieve good agreement in recovering the wave elevation, albeit with the necessity of using refined grids in the regions of greatest surface variation. Refining a grid where needed is far much easier than finding a better basis of functions, that is, a feature of considerable practical interest.

The method presented in this study consists of the first boundary integral approach to solve this nonlinear recovery problem. We expect this work to pave the way for addressing even more challenging configurations, including the extension to three-dimensional settings, which will inevitably involve Green functions in the integral kernels.

The paper is organised as follow. The mathematical model and relations of interest are introduced in § 2, with an Eulerian description of motion. In order to handle overhanging waves, their Lagrangian counterparts are introduced in § 4. In § 3 we derive an equation for the free surface, allowing us to compute reference solutions for testing our recovery procedure. This procedure is described in the subsequent § 5. Numerical implementation and examples are provided in § 6. Finally, in § 7 we discuss our general results, as well as possible future extensions and implications of this work.

2. Mathematical settings

We give here the classical formulation of our water-wave problem in Eulerian variables. Physical assumptions and notations being identical to that of Clamond et al. (2023), interested readers should refer to this paper for further details.

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. 2023). 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. 2023):

(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. 2023). 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. (2023) for more details).

Following Clamond & Constantin (2013), 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 2006; Constantin & Strauss 2010).

As for irrotational waves, it is useful (Clamond 2013, 2018; Clamond et al. 2023) 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 1953), 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 1990)

(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.

3. Lagrangian description

This section focuses on addressing the limitation of the Eulerian framework that hinders the computation of overhanging waves, which are characterised by multi-valued surfaces. While one option to address this challenge involves employing an arclength formulation, as elaborated by Vanden-Broeck (1994), we have chosen to adopt a Lagrangian formalism in this study. One benefit of the Lagrangian approach is its inherent capability to aggregate collocation points near the wave crest, where they are most crucial (Da Silva & Peregrine 1988).

Our approach is similar to that used by Da Silva & Peregrine (1988), although we modify the integral kernels in terms of polylogarithms and we integrate the Cauchy integral formula to remove the strong singularity from the kernels following Clamond (2018). The resulting analytical expression is characterised by a weak logarithmic singularity, and it is suitable for calculating various types of waves, including solitary and periodic waves, overhanging or not.

Considering the (rotational) complex velocity $w=u-\mathrm {i} v$ evaluated at the surface, let us introduce the holomorphic function $\log {(-{w}/\sqrt {gd})}={q}-\mathrm {i}{\theta }$ and define accordingly

(3.1)$$\begin{gather} {q}_{s} \stackrel{\text{def}}= \operatorname{Re}{\{\log{(-{w}_{s}/\sqrt{gd})}\}} = \tfrac{1}{2} \ln{ [ ({u}_{s}^2+{v}_{s}^2)/\sqrt{gd} ]}, \end{gather}$$

(3.2)$$\begin{gather}{\theta}_{s} \stackrel{\text{def}}={-}\operatorname{Im}{\{\log{(-{w}_{s}/\sqrt{gd})}\}} ={-}\operatorname{atan2}{({v}_{s},-{u}_{s})}, \end{gather}$$

where we notably used the Bernoulli principle at the surface and where $\operatorname {atan2}$ denotes the four-quadrant inverse tangent.

Exploiting the impermeability and isobarity of the free surface, one gets

(3.3)\begin{equation} \tan{{\theta}_{s}} = \eta_x = \frac{{v}_{s}}{{u}_{s}} = \frac{\sigma{v}_{s}}{\sqrt{{B}_{s} - 2g\eta - {v}_{s}^2}} = \frac{\sqrt{{B}_{s} - 2g\eta - {u}_{s}^2}}{\sigma{u}_{s}}, \end{equation}

where $\sigma =\pm 1$ is introduced for conveniently choosing the direction of the wave propagation (see the discussion at the beginning of § 4.1 below). Hence, extracting ${u}_{s}$ and ${v}_{s}$ from the latter relations, we have

(3.4a,b)\begin{equation} {u}_{s} = \sigma\cos({{\theta}_{s}}) \sqrt{{B}_{s} - 2g\eta},\quad {v}_{s} = \sigma\sin({{\theta}_{s}}) \sqrt{{B}_{s} - 2g\eta}. \end{equation}

We now consider the Lagrangian description of motion, with $t$ denoting the time, thus, ${u}_{s}=\mathrm {d}\kern 0.06em x/\mathrm {d} t$ and ${v}_{s}=\mathrm {d}\eta /\mathrm {d} t$ ($\mathrm {d}/ \mathrm {d} t$ is the temporal derivative following the motion). From the second expression in (3.4a,b), we deduce that

(3.5)\begin{equation} \frac{\textrm{d}}{\textrm{d} t} \sqrt{{B}_{s} - 2g\eta} ={-} \sigma g \sin({{\theta}_{s}}), \end{equation}

and hence, considering a crest at $t=0$ (where $\eta (0)=a$ is the wave amplitude), we have

(3.6a,b)\begin{equation} \sqrt{{B}_{s} - 2g\eta} = \mu - g \sigma \int_{0}^{t} \sin({{\theta}_{s}'})\, \textrm{d} t',\quad \mu \stackrel{\text{def}}= \sqrt{{B}_{s} - 2ga}, \end{equation}

where ${\theta }_{s}'\stackrel {\small {\rm def}}={\theta }_{s}(t')$. Therefore, all quantities at the free surface can be expressed in terms of the surface angle ${\theta }_{s}$, using $t$ as an independent variable, e.g.

(3.7)$$\begin{gather} {z}_{s}(t) = \sigma \int_{0}^{t} \Bigg[ \mu - g \sigma \int_{0}^{t'} \sin({{\theta}_{s}''})\, \textrm{d} t'' \Bigg] \exp{(\textrm{i}{\theta}_{s}')}\,\textrm{d} t' + \textrm{i} a , \end{gather}$$

(3.8)$$\begin{gather}{w}_{s}(t) = \frac{\textrm{d}{\bar{z}}_{s}}{\textrm{d} t} = \sigma \left[ \mu - g \sigma \int_{0}^{t} \sin({\theta}_{s}')\,\textrm{d} t' \right] \exp{(-\textrm{i}{\theta}_{s})}, \end{gather}$$

(3.9)$$\begin{gather}\eta(t) = \operatorname{Im}{{z}_{s}} = \frac{{B}_{s}}{2g} - \frac{1}{2g} \left[ \mu - g \sigma \int_{0}^{t} \sin({\theta}_{s}')\,\textrm{d} t' \right]^2 . \end{gather}$$

In the Eulerian description of motion, the wave period is constant and depends on the reference frame (moving at a constant phase speed) where one observes the fluid. On the other hand, in the Lagrangian context, the period $T_L$ of the free surface differs from the Eulerian period due to the Stokes drift (Longuet-Higgins 1986). The Lagrangian period being such that $\eta (t+T_L)=\eta (t)$, we exploit expression (3.9) which, after some elementary algebra and simplifications, yields

(3.10)\begin{equation} \left[ 2\mu - g \sigma \int_{t}^{t+T_L} \sin({\theta}_{s}')\, \textrm{d} t' - 2g\sigma \int_{0}^{t} \sin({\theta}_{s}')\, \textrm{d} t'\right] \int_{t}^{t+T_L} \sin({\theta}_{s}')\, \textrm{d} t' = 0, \end{equation}

that is necessarily satisfied for all times if and only if

(3.11)\begin{equation} \int_{0}^{T_L} \sin({\theta}_{s}')\, \textrm{d} t' = 0, \end{equation}

thus defining the Lagrangian period $T_L$.

4. Equations for the free surface

Since we do not have access to measurements of the bottom pressure for nonlinear waves with constant vorticity, we must generate these data from numerical solutions of the exact equations. Thus, this section aims to derive a comprehensive formulation for the computation of surface waves using a boundary integral method. From these solutions, the bottom pressure is subsequently obtained and used as input for our surface recovery procedure.

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 2018)

(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 (2018), 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 (1988) and by Vanden-Broeck (1994) – thus allowing us to further reduce the numerical cost.

5. Surface recovery from bottom pressure

Building upon the last three sections, we propose a nonlinear integral equation to recover the surface elevation in terms of a given measure of pressure ${p}_{b}(x)$ at the seabed. This ‘measurement’ is given here numerically by solving expression (4.10) and, subsequently, computing the pressure from the surface profile (exploiting the boundary integrals and the Bernoulli equation).

The principal benefit in establishing an integral formulation lies in the fact that no specific eigenfunctions are needed to fit the pressure data. In fact, our derived equation remains applicable regardless of the nature of the wave under consideration, be it periodic or not, overhanging or not. We re-emphasise here that, although the equations consider $(2{\rm \pi} /k)$-periodic waves, their aperiodic counterparts are easily obtained letting $k\to 0^+$.

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. Numerical illustrations

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 2013; Clamond et al. 2023). 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. (2023). 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 (1988), 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 2001), 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 1971). 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 (1988). 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.

7. Discussion

This work presents a novel and comprehensive boundary integral method for recovering surface water waves from bottom pressure measurements. Despite the inherent complexity of this inverse problem, we successfully formulate the relatively simple expression (5.9) for surface recovery, enabling the computation of a wide range of rotational steady waves. A significant advantage of this approach lies in the integral formulation, which eliminates the need for arbitrarily selecting a basis of functions to fit the pressure data, as done previously (Clamond 2013; Clamond & Constantin 2013; Clamond et al. 2023).

To demonstrate the robustness and efficiency of our method, we showcase two challenging examples: an overturning wave with a large amplitude and a solitary wave. In both cases, we accurately recovered the surface profile and the hydrodynamic parameters with good agreement. Although it might be possible to adapt our numerical procedure to compute extreme waves (with angular surface) with (or without) overhanging profiles, this task is left for future investigations. In fact, our main goal here is a proof of concept and to provide clear evidence on the effectiveness of this new formulation. It is worth noting that the existence of solitary wave solutions for this inverse problem was established by Henry (2013) for steady waves with analytic vorticity distributions. It still remains an open question whether this statement holds or not for periodic waves.

In conclusion, this paper, along with the proposed boundary integral formulation, represents a significant milestone in solving the surface wave recovery problem, providing a solid foundation for future extensions, such as its potential application to three-dimensional configurations. Indeed, in three-dimensional holomorphic functions cannot be employed but integral representations via Green functions remain, so an efficient fully nonlinear surface recovery is conceivable.