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 Reference Tsai and Tsai2009). 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 Reference Lin and Yang2020). 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 Reference Bergan, Tørum and Traetteberg1969; Lee & Wang Reference Lee and Wang1985), it was only recently that nonlinear waves began to be addressed (Oliveras et al. Reference Oliveras, Vasan, Deconinck and Henderson2012). Some works, such as Constantin (Reference Constantin2012), 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 Reference Clamond2013; Clamond & Constantin Reference Clamond and Constantin2013). These studies demonstrated the convergence of the reconstruction process and the ability to recover waves of maximum amplitude (Clamond & Henry Reference Clamond and Henry2020). Furthermore, the method was adapted to handle cases involving linear shear currents (Clamond, Labarbe & Henry Reference Clamond, Labarbe and Henry2023), 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 Reference Wahlen2009). Note that the existence of nonlinear water waves with constant vorticity and overturning wave profiles has been demonstrated (Constantin & Varvaruca Reference Constantin and Varvaruca2011; Constantin, Strauss & Varvaruca Reference Constantin, Strauss and Varvaruca2016). Note also that, for infinitesimal waves, exact recovery formulae were obtained by Henry & Thomas (Reference Henry and Thomas2018) for steady waves with general vorticity.
The recovery method studied in Clamond (Reference Clamond2013), Clamond & Constantin (Reference Clamond and Constantin2013), Clamond & Henry (Reference Clamond and Henry2020) and Clamond et al. (Reference Clamond, Labarbe and Henry2023) 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 Reference Clamond2013; Clamond & Constantin Reference Clamond and Constantin2013). 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 (Reference Da Silva and Peregrine1988) 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 Reference Vanden-Broeck1994), 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. (Reference Clamond, Labarbe and Henry2023), 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
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.
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
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
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
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):
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,
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
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
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)$,
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.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
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
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)
Alternatively, using the method of images, along with the identity (A10), the Cauchy integral can be written as
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
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
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 (Reference Vanden-Broeck1994), 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 Reference Da Silva and Peregrine1988).
Our approach is similar to that used by Da Silva & Peregrine (Reference Da Silva and Peregrine1988), 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 (Reference Clamond2018). 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
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
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
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
and hence, considering a crest at $t=0$ (where $\eta (0)=a$ is the wave amplitude), we have
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.
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 Reference Longuet-Higgins1986). 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
that is necessarily satisfied for all times if and only if
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
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),
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
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
where we exploited the property (Clamond Reference Clamond2018)
Similarly, we have
Finally, after integrating the whole expression (4.3) and retaining the imaginary part only, we obtain the equation for the free surface
where $K$ is an integration constant obtained enforcing the mean-level condition (2.1), i.e.
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
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
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.
where we introduced
with $\tau \stackrel{\text{def}}{=} 2{\rm \pi} /T_L$. Considering $t\to t'$ in both integrands in (4.11), we have
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. 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
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
Before continuing, let the compact notations for the polylogarithmic kernels be
We now evaluate expression (5.2) at the free surface, reducing it to
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,
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.
After some algebraic manipulations, it further reduces to
Substituting (5.7) into (5.5), the imaginary part yields the Eulerian integral formulation for the surface recovery
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:
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 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
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}$.
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.
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
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.
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 Reference Clamond2013; Clamond & Constantin Reference Clamond and Constantin2013; Clamond et al. Reference Clamond, Labarbe and Henry2023).
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 (Reference Henry2013) 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.
Acknowledgements
The authors would like to thank anonymous referees for comments that helped to improve the paper readability.
Funding
J.L. has been supported by the French government, through the $\mbox {UCA}^{JEDI}$ Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Logarithms and polylogarithms
The function $\ln (x)$ denotes the natural logarithm (Napierian logarithm) of a real positive variable $x$ ($x\in \mathbb {R}^+$), and $\log (z)$ denotes the principal logarithm of a complex variable $z\in \mathbb {C}$, i.e.
This definition requires that the argument of any complex number lies in $]-{\rm \pi} ;{\rm \pi} ]$. It implies, in particular, that $\arg (z^{-1})=-\arg (z)$ if $z\not \in \mathbb {R}^-$ and that $\arg (z^{-1})=\arg (z)={\rm \pi}$ if $z\in \mathbb {R}^-$. We have the special relations
where $\lfloor {\cdots }\rfloor$ is the rounding toward $-\infty$.
The polylogarithms can be defined, for $|z|<1$ and $\nu \in \mathbb {C}$, by
and for all complex $z$ by analytic continuation (Wood Reference Wood1992). With the above definition of the complex logarithm, we have the special inversion formulae
The polylogarithms for $\nu \in \mathbb {N}^*$ are single-valued functions in the cut plane $z\in \mathbb {C}\backslash [1;+\infty [$ and the inversion formula can be used to extend their definition for $z\in [1;+\infty [$. Note that the inversion formula depends on the definition of the principal logarithm that is not unique, thus, several variants can be found in the literature.
We finally note other relations useful in this paper,
where $[{\cdots }]$ denotes the rounding toward zero, the last relation being valid for $-2{\rm \pi} <\mathrm {Re}\{z\}<2{\rm \pi}$.
Appendix B. Singular Leibniz integral rule
Let $K(x,y)$ be a function regular everywhere for $y\in [a,b]$, except perhaps at $y=y_0\in ]a,b[$ ($y_0$ generally depending on $x$, as well as $a$ and $b$) where $K$ may be singular, its finite integral being taken in the sense of Cauchy's principal value, i.e.
exists. The first derivative of $J$ is thus
For instance, $a$ and $b$ being constant and for a sufficiently well-behaving function $\varphi$, we have
where the limit is zero if $\varphi$ is Hölder continuous (sufficient but unnecessary condition). This formula is a consequence of (B2) but also of the choice of the antiderivative of $1/(x-y)$. Indeed, one can also write
where we used the relations $\log (\epsilon )=\ln (\epsilon )$ and $\log (-\epsilon )= \ln (\epsilon )+\mathrm {i}{\rm \pi}$, since $\epsilon \in \mathbb {R}^+$. The relation (B3) is more convenient when dealing only with real variables, while (B4) is more suitable for complex formulations. The reason behind the latter argument is because $\log (x)$ can be continued analytically in the complex plane, which is not the case with $\ln |x|$.