# Optimal reconstruction of water-waves from noisy pressure measurements at
the seabed

Joris Labarbe  Université Côte d’Azur, CNRS UMR 7351, Laboratoire J. A.
Dieudonné, Parc Valrose, 06108 Nice cedex 2, France  joris.LABARBE@univ-
cotedazur.fr ,  Alexandre Vieira  Université Côte d’Azur, CNRS UMR 7351,
Laboratoire J. A. Dieudonné, Parc Valrose, 06108 Nice cedex 2, France
<EMAIL_ADDRESS>and  Didier Clamond  Université Côte
d’Azur, CNRS UMR 7351, Laboratoire J. A. Dieudonné, Parc Valrose, 06108 Nice
cedex 2, France<EMAIL_ADDRESS>
###### Abstract.

We consider the problem of recovering the surface wave profile from noisy
bottom pressure measurements with (a priori unknown) arbitrary pressure at the
surface. Without noise, the direct approach developed in [clamond2023steady]
provides an effective way to recover the sea surface. However, the assumption
of analyticity for the measurement renders this method inefficient in the
presence of noise. Therefore, we introduce an optimisation procedure based on
the minimisation of a distance between a recovered bottom pressure and its
measurement. Such method proves to be well-designed to handle perturbed
signals. We illustrate the effectiveness of this approach in the recovery of
gravity-capillary waves from unfiltered noisy data.

######  Mathematics Subject Classification:

##  1\. Introduction

Monitoring the surface of the oceans is a timely concern for climatic and
environmental considerations, especially in coastal regions where large waves
represent a risk when approaching the shoreline. Currently, an efficient way
of monitoring the wave motion is lacking, as most current techniques display a
poor accuracy in their measurements. For this reason, in the second half of
the last century, some scientists emitted the idea of sending probes at the
seabed to measure the pressure and reconstruct the surface profile from these
data [BTT69]. This procedure allows to determine the shape of water-waves
without intrusive observations, although it necessitates to solve an inverse
problem within an unknown domain. Such problem is notoriously challenging,
notably because of its ill-posedness (any disturbance grows exponentially from
the bottom to the free surface).

Historically, the first attempt to solve this problem relied on hydrostatic
theory and assumed the surface elevation as being proportional to the weight
of the water column. Later, from a linear approximation of the problem, an
explicit expression for the surface was given in terms of a certain pressure
transfer function [LW85]. Noticeably, a recent work derived an exact formula
for steady rotational waves using linear theory [HT16]. However, both these
assumptions were not suited to describe large waves or of non-trivial profile
(i.e. not sinusoidal). It was only recently that the first results on the full
nonlinear problem were obtained [O12]. This study resulted in a nonlocal
integro-differential equation for the surface elevation that was highly
difficult to solve in practice without considering asymptotic limits. Slightly
after this work, an efficient procedure for the recovery of nonlinear steady
water-waves was formulated [C13, CC13]. This approach relies on the assumption
that bottom pressure can be analytically continued as a holomorphic function,
allowing to express the surface elevation as a solution of an algebraic
equation. Numerical investigation of this method was done for extreme waves
[CC13] and it was rigorously proved to converge in the case of steady surface
gravity waves from a mere fixed-point algorithm [CH20]. Subsequently, this
theory was extended to handle more complex configurations, such as the
presence of a linear shear current [CLH23], the possibility of overhanging
waves [LC23] or the influence of arbitrary pressure at the surface
[clamond2023steady]. We note the proof for the existence of nonlinear water-
waves with constant vorticity and overturning wave profiles [CSV11, CSV16].

The method originally developed in [CC13] consistently proved its
effectiveness in terms of accuracy, computation time and accessibility, even
considering technical difficulties (stagnation points, limiting waves,
capillary effects, etc.). However, the formulation inherently depends on the
analyticity of the measurement, which is not the case in real-life
observations where the presence of external noise is inevitable. Moreover,
since we are performing most of the calculations in the complex plane, then it
prohibits us to directly extend this approach to the three-dimensional
configuration. For all these reasons, we introduce a new and innovative
formulation, based on an optimisation method, to solve this constrained
inverse problem with noise in the measurements.

To the best of authors’ knowledge, surface wave recovery within an
optimisation framework has never been considered in the literature. As it
turns out, one can relate this reconstruction problem to an optimal control
approach, where the goal is to find a set of parameters or functions
minimising a given cost functional under some dynamical constraints. (We refer
interested readers to [cossu2014introduction] for an introduction on this
subject.) In the current inverse problem, the undetermined function is the
surface profile and the cost functional will be the distance to the
observation. In this sense, this problem is close to some shape optimisation
or topology optimisation (TO) configurations. These problems consist in
finding the optimal shape of an obstacle immersed in a moving fluid (see
[alexandersen2020review, dbouk2017review] and references therein for an
overview on the numerical resolution of TO problems applied to several
physical settings). Although it yields satisfactory results in general, the
main drawback of the shape optimisation problem comes from its numerical
complexity. This is because the underlying dynamical constraints are usually
discretised using a Finite Element method, which relies on a prescribed mesh
adapted to the shape of the obstacle. In this setting, the shape of the
immersed obstacle is updated at each iteration, leading to an expensive
recalculation of the mesh [allaire2014shape]. Thus, all techniques avoiding
this costly procedure necessarily enhance the formulation.

The method we consider in this work relies on a predefined discretisation of
some analytical expressions we derived in an earlier work, from a boundary
integral method [LC23]. Combining this analytic preprocessing within the
optimisation framework has never been explored before and appears to be a
suitable and innovative approach for the surface reconstruction problem.
Moreover, the algorithm we develop is easily extendable to a broad range of
situations (e.g. adding further physical assumptions or noise in the
measurements) and yields an excellent agreement with the sought solution
(comparable to the results obtained in [clamond2023steady] using the
holomorphic extension).

This article is articulated around the presentation, application and
comparison of two different methods to reconstruct surface waves when
considering noise or not. Section 2 introduces the equations of motion and
boundary conditions for the system of interest. Then, sections 3 and 4
describe in details the direct and optimisation approaches, respectively, and
how to implement them numerically in an efficient way. A numerical example is
treated in section 5 to explore and compare thoroughly these methods. Finally,
we discuss about the results we obtained and on the eventual improvements and
future extensions in section 6.

##  2\. Problem statement

We investigate the motion of a steady wave travelling at a constant phase
speed $c$ within a fluid layer with a free surface. The bottom of the layer
being flat, we fix ourselves in the Galilean frame of reference moving with
the wave, so that motion appears stationary for the observer. The fluid is
considered inviscid, irrotational and incompressible; we assume the surface
and seabed to be impermeable. These boundaries are located at $y=\eta(x)$ for
the surface and $y=-0pt$ for the seabed, with $y$ being the vertical
coordinate and $0pt$ a constant depth. We assume the streamwise direction to
have an infinite extension. Therefore, the domain of definition is given by
$\Omega=\\{(x,y):x\in\mathds{R},-0pt\leqslant y\leqslant\eta(x)\\}$. When the
motion is periodic, we introduce the wave period $L=2{\pi}/k$ (the limit $k\to
0^{+}$ corresponding to a solitary wave) where $k$ is a spatial wavenumber.
Thus, we define an Eulerian averaging operator to fix the mean water level at
$y=0$, i.e.

(1)
$\left<\eta\right>\vcentcolon=\frac{k}{2{\pi}}\int_{\mathcal{P}}\eta(x)\mathrm{d}x=0,$

where the path of integration is $\mathcal{P}=[-{\pi}/k,{\pi}/k]$.

In the following, we use the subscripts ‘s’ and ‘b’ to denote the restriction
of fields at the surface and the bottom, respectively. Alternatively, the
subscript $x$ is used to express horizontal derivatives (e.g.
$\eta_{x}\vcentcolon=\mathrm{d}\eta/\mathrm{d}x)$. We consider the pressure at
the surface as an undetermined function of the $x$-coordinate, entering the
equations through the dynamical condition (supplemented with kinematic
condition). The velocity field $\bm{u}=(u,v)$ is governed by the steady Euler
equations (expressing the conservation of mass and linear momentum) with
associated boundary conditions. It yields

(2a) $\displaystyle\bm{\nabla}\bm{\cdot}\bm{u}$ $\displaystyle=0,\quad$
$\displaystyle\textrm{in}\quad\Omega$ (2b)
$\displaystyle\bm{u}\bm{\cdot}\bm{\nabla}\bm{u}+\bm{\nabla}p+\bm{g}$
$\displaystyle=0,\quad$ $\displaystyle\textrm{in}\quad\Omega$ (2c)
$\displaystyle\bm{u}\bm{\cdot}\bm{n}-\bm{u}\bm{\cdot}\bm{\nabla}\eta$
$\displaystyle=0\quad$ $\displaystyle\textrm{at}\quad y=\eta(x),$ (2d)
$\displaystyle p$ $\displaystyle={p}_{\mathrm{s}}(x)\quad$
$\displaystyle\textrm{at}\quad y=\eta(x),$ (2e)
$\displaystyle\bm{u}\bm{\cdot}\bm{n}$ $\displaystyle=0\quad$
$\displaystyle\textrm{at}\quad y=-d,$

where $p$ is a relative pressure with zero mean-value at the surface (i.e.
$\langle{p}_{\mathrm{s}}\rangle=0$) and scaled with the constant background
density $\rho_{0}$. The unit normal vector $\bm{n}$ is directed outward and
$\bm{g}=(0,g)$ represents the restoring gravity force, with $g$ the
acceleration due to gravity (acting downwards). By integrating equation (2b),
we recover the well-known Bernoulli principle (expressing the conservation of
mechanical energy), leading (in the irrotational case) to

(3) $2(p+gy)+u^{2}+v^{2}=B,$

where $B$ is the Bernoulli constant. From this principle and using the mean
operator (1), we establish some important relations, notably that
$\langle{u}_{\mathrm{s}}^{2}+{v}_{\mathrm{s}}^{2}\rangle=B$ and
$\langle{p}_{\mathrm{b}}\rangle=g0pt$. Finally, we introduce the complex
velocity $w(z)\vcentcolon=u(x,y)-\mathrm{i}v(x,y)$ that is holomorphic in the
$z$-coordinate (with $z\vcentcolon=x+\mathrm{i}y$). We note that the free
surface path is given by the curvilinear abscissa
${z}_{\mathrm{s}}=x+\mathrm{i}\eta$ whereas the solid boundary is described by
${z}_{\mathrm{b}}=x-\mathrm{i}0pt$.

The principal objective of this work is to recover the surface wave profile
$\eta(x)$, the Bernoulli constant $B$ and the surface pressure
${p}_{\mathrm{s}}(x)$ from a given pressure observation
${p}_{\mathrm{b}}^{\mathrm{obs}}(x)$ at the seabed. We assume either the wave
period $L$ or the wavenumber $k$ to be a known parameter. This inverse problem
is notoriously difficult to solve due to the strong nonlinearities and its
ill-posed nature (i.e. disturbances are exponentially growing to the surface).
Our first approach, developed in a series of articles [C13, CC13, CH20,
clamond2023steady, CLH23, LC23], is a direct approach set in the physical
space and where analytical expressions (integro-differential in general) are
given in terms of the unknown surface wave profile and surface pressure.

##  3\. Numerical procedure: direct approach

###  3.1. Cauchy integral formula

Computing steady surface waves from a boundary integral formula is an idea
initiated in a seminal work by Da Silva and Peregrine [DSP88]. Hereafter, we
present a similar procedure to recover an implicit expression for the surface
elevation using complex analysis [LC23]. We start by giving the Cauchy
integral formula, written for a holomorphic function $\Upxi(z)$, in its
classical form

(4)
$\mathrm{i}\vartheta\Upxi(z)=\text{P.V.}\oint\frac{\Upxi(z^{\prime})}{z^{\prime}-z}\mathrm{d}z^{\prime}=\int_{-\infty}^{\infty}\frac{{\Upxi}_{\mathrm{b}}^{\prime}\mathrm{d}x^{\prime}}{{z}_{\mathrm{b}}^{\prime}-z}-\int_{-\infty}^{\infty}\frac{\left(1+\mathrm{i}\eta_{x}^{\prime}\right){\Upxi}_{\mathrm{s}}^{\prime}\mathrm{d}x^{\prime}}{{z}_{\mathrm{s}}^{\prime}-z},$

with primes denoting the dependence on the dummy variable (e.g.
$\Upxi^{\prime}\vcentcolon=\Upxi(x^{\prime})$). The argument $\vartheta$ is
substituted by $\\{2\pi,0,\pi\\}$ when the coordinate $z$ lies respectively
inside, outside and at the smooth boundary of the domain $\Omega$.

Whenever the holomorphic function is purely real at the bottom boundary (i.e.
$\operatorname{Im}\\{{\Upxi}_{\mathrm{b}}\\}=0$), we may use the method of
images (also called Schwarz reflection principle [MorseFeshbach1953]) on (4)
to obtain

(5)
$\displaystyle\vartheta\Upxi(z)=k\int_{\mathcal{C}}\mbox{L{\scriptsize{}}}_{0}\\{\mathrm{e}^{\mathrm{i}k({z}_{\mathrm{s}}^{\prime}-z)}\\}{\Upxi}_{\mathrm{s}}^{\prime}\mathrm{d}{z}_{\mathrm{s}}^{\prime}+k\int_{\mathcal{C}}\mbox{L{\scriptsize{}}}_{0}\\{\mathrm{e}^{\mathrm{i}k(z-{\tilde{z}}_{\mathrm{s}}^{\prime}+2\mathrm{i}0pt)}\\}{\widetilde{\Upxi}}_{\mathrm{s}}^{\prime}\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}+\pi\left<{\Upxi}_{\mathrm{s}}\frac{\mathrm{d}{z}_{\mathrm{s}}}{\mathrm{d}x}+{\widetilde{\Upxi}}_{\mathrm{s}}\frac{\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}}{\mathrm{d}x}\right>,$

where $\mbox{L{\scriptsize{}}}_{\nu}$ is the $\nu$th polylogarithm [LC23] and
$\widetilde{z}$ the complex conjugate of $z$. The integration path is given
(in the complex plane) by $\mathcal{C}=\\{{z}_{\mathrm{s}}(x)\ :\
x\in[-\pi/k,\pi/k]\\}\subset\mathbb{C}$.

Equation (5) possesses a strong (polar) singularity causing technical
difficulties when computing the integral numerically. Thus, we use the
integrated form of the polylogarithm to replace the polar term with a
logarithmic singularity [LC23]. We obtain

(6)
$\displaystyle\vartheta\Upxi(z)=\mathrm{i}\\!\int_{\mathcal{C}}\partial_{z}\mbox{L{\scriptsize{}}}_{1}\\{\mathrm{e}^{\mathrm{i}k({z}_{\mathrm{s}}^{\prime}-z)}\\}{\Upxi}_{\mathrm{s}}^{\prime}\mathrm{d}{z}_{\mathrm{s}}^{\prime}+\mathrm{i}\\!\int_{\mathcal{C}}\partial_{z}\mbox{L{\scriptsize{}}}_{1}\\{\mathrm{e}^{\mathrm{i}k(z-{\tilde{z}}_{\mathrm{s}}^{\prime}+2\mathrm{i}0pt)}\\}{\widetilde{\Upxi}}_{\mathrm{s}}^{\prime}\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}+\pi\left<{\Upxi}_{\mathrm{s}}\frac{\mathrm{d}{z}_{\mathrm{s}}}{\mathrm{d}x}+{\widetilde{\Upxi}}_{\mathrm{s}}\frac{\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}}{\mathrm{d}x}\right>,$

Computing the latter expression at the free surface reduces to

(7)
$\displaystyle{\Upxi}_{\mathrm{s}}\mathrm{d}{z}_{\mathrm{s}}=\frac{\mathrm{i}}{2\pi}\mathrm{d}\\!\left[\int_{\mathcal{C}}\mbox{L{\scriptsize{}}}_{1}\\{\mathrm{e}^{\mathrm{i}k({z}_{\mathrm{s}}^{\prime}-z)}\\}{\Upxi}_{\mathrm{s}}^{\prime}\mathrm{d}{z}_{\mathrm{s}}^{\prime}+\int_{\mathcal{C}}\mbox{L{\scriptsize{}}}_{1}\\{\mathrm{e}^{\mathrm{i}k(z-{\tilde{z}}_{\mathrm{s}}^{\prime}+2\mathrm{i}0pt)}\\}{\widetilde{\Upxi}}_{\mathrm{s}}^{\prime}\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}\right]+\frac{\mathrm{d}{z}_{\mathrm{s}}}{2}\left<{\Upxi}_{\mathrm{s}}\frac{\mathrm{d}{z}_{\mathrm{s}}}{\mathrm{d}x}+{\widetilde{\Upxi}}_{\mathrm{s}}\frac{\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}}{\mathrm{d}x}\right>.$

###  3.2. Expressions for computing steady surface waves and bottom pressure

From Bernoulli’s principle (3), the complex velocity at the surface
${w}_{\mathrm{s}}$ is given explicitly by

(8)
${w}_{\mathrm{s}}=\sigma(\mathrm{d}\widetilde{{z}_{\mathrm{s}}}/\mathrm{d}x)\sqrt{(B-2{p}_{\mathrm{s}}-2g\eta)/|\mathrm{d}{z}_{\mathrm{s}}/\mathrm{d}x|^{2}},$

where $\sigma=\mp 1$ denotes travelling waves propagating along or against the
current, respectively.

Let us consider the holomorphic function $\Upxi=w+c$ ($c$ being an arbitrary
definition of the phase speed). The left-hand side of (7) follows directly
from (8) as

(9) ${\Upxi}_{\mathrm{s}}\mathrm{d}{z}_{\mathrm{s}}=\left(\omega
h+c\right)\mathrm{d}{z}_{\mathrm{s}}+\sigma\sqrt{(B-2{p}_{\mathrm{s}}-2g\eta)/|\mathrm{d}{z}_{\mathrm{s}}/\mathrm{d}x|^{2}}\mathrm{d}x,$

where $h\vcentcolon=\eta+0pt$ is the wave height. We note that the radicand is
purely real since $B\geqslant 2||{p}_{\mathrm{s}}+g\eta||_{\infty}$ for all
waves (according to Bernoulli’s relation).

After integration of expression (7) – retaining the imaginary part only – we
obtain an equation for the computation of the free surface

(10) $\frac{\omega\/\eta^{2}}{2}-\eta\left<{u}_{\mathrm{s}}\left|\\\
ud{z}_{\mathrm{s}}/\mathrm{d}x\right|^{2}\right>-\frac{\sigma}{2\pi}\int_{\mathcal{P}}\operatorname{Re}\\{\mathcal{L}_{1}\\}\sqrt{(B-2{p}_{\mathrm{s}}^{\prime}-2g\eta^{\prime})/|\mathrm{d}{z}_{\mathrm{s}}^{\prime}/\mathrm{d}x^{\prime}|^{2}}\mathrm{d}\/x^{\prime}-K=0,$

where $K$ is a constant of integration, recovered by enforcing the mean-level
condition [LC23]. For the sake of brevity, we also introduced the following
notation

$\mathcal{L}_{\nu}\vcentcolon=\mbox{L{\scriptsize{}}}_{\nu}\\!\big{[}\mathrm{e}^{\mathrm{i}k({z}_{\mathrm{s}}^{\prime}-{z}_{\mathrm{s}})}\big{]}-\mbox{L{\scriptsize{}}}_{\nu}\\!\big{[}\mathrm{e}^{\mathrm{i}k({z}_{\mathrm{s}}-{\tilde{z}}_{\mathrm{s}}^{\prime}+2\mathrm{i}0pt)}\big{]}.$

As one can notice, expression (10) is not dependent on the choice of phase
speed $c$, as expected from the Galilean invariance of the problem.

Equation (5) and the Schwarz reflection principle also let us find an
expression for the bottom velocity. Consider $\Upxi=w$ as the holomorphic
function in (5) and evaluate the expression at the bottom (with now
$\vartheta=2{\pi}$). Since ${u}_{\mathrm{b}}$ is a real function of the
horizontal coordinate, we use the method of images to recover its explicit
form, which leads to the formula

(11)
${u}_{\mathrm{b}}(x)=\frac{\mathrm{i}k}{4\pi}\left[\int_{\mathcal{C}}\cot\\!\left(k\frac{{z}_{\mathrm{s}}^{\prime}-{z}_{\mathrm{b}}}{2}\right){w}_{\mathrm{s}}^{\prime}\mathrm{d}{z}_{\mathrm{s}}-\int_{\mathcal{C}}\cot\\!\left(k\frac{{\tilde{z}}_{\mathrm{s}}^{\prime}-{z}_{\mathrm{b}}-2\mathrm{i}0pt}{2}\right){\widetilde{w}}_{\mathrm{s}}^{\prime}\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}\right].$

This bottom velocity ${u}_{\mathrm{b}}$ is used in Stokes’ first definition of
the phase speed

(12) $c_{1}=-\left<{u}_{\mathrm{b}}\right>.$

###  3.3. Holomorphic functions and recovery formula

As done numerous time in previous articles [C13, CC13, CH20,
clamond2023steady, CLH23, LC23], we introduce a holomorphic complex pressure
function $\mathfrak{P}(z)$ by analytic continuation in the domain $\Omega$.
Using Bernoulli’s principle written for the complex velocity, it yields an
expression for the complex pressure as

(13) $\mathfrak{P}(z)\vcentcolon=g0pt+\frac{B-w^{2}}{2}.$

Similarly, we introduce the antiderivative of $\mathfrak{P}$ defined by

(14)
$\displaystyle\mathfrak{Q}(z)\vcentcolon=\int_{z_{0}}^{z}\left[\mathfrak{P}(z^{\prime})-g0pt\right]\mathrm{d}z^{\prime}=\frac{1}{2}\int_{z_{0}}^{z}\left[B-w(z^{\prime})^{2}\right]\mathrm{d}z^{\prime},$

where $z_{0}$ is an arbitrary constant.

In practice, these holomorphic functions are easily recovered from the given
pressure measurement by fitting the data over a suitable eigenbasis (one
leading to the most accuracy with the lowest order of quadrature). Assuming a
periodic motion, we fit these data with a Fourier polynomial basis before
performing an analytic continuation of the pressure field. Evaluated at the
surface, it yields

(15) $\mathfrak{P}(z)={p}_{\mathrm{b}}(z+\mathrm{i}0pt)\approx
g0pt+\sum_{|n|>0}^{N}\mathfrak{p}_{n}\mathrm{e}^{\mathrm{i}nk(z+\mathrm{i}d)}.$

From this definition, we simply integrate once and recover the second
holomorphic function as

(16)
${\mathfrak{Q}}_{\mathrm{s}}(x)=\int_{0}^{x}\left[{\mathfrak{P}}_{\mathrm{s}}(x^{\prime})-g0pt\right]\mathrm{d}x^{\prime}\approx\sum_{|n|>0}^{N}\frac{\mathrm{i}\mathfrak{p}_{n}}{nk}\frac{\mathrm{e}^{-nka}-\mathrm{e}^{\mathrm{i}nk(x+\mathrm{i}\eta)}}{\mathrm{e}^{nkd}},$

where $a$ is the amplitude of the wave at the crest (arbitrarily located at
$x=0$).

In the context of the direct approach, we exploit a relation derived in a
recent work as our general recovery formula [clamond2023steady]. This
expression (not involving any differential terms of $\eta$) is given by

(17)
$\operatorname{Re}\\!\left\\{{\mathfrak{Q}}_{\mathrm{s}}\right\\}=\frac{1}{2}\int_{0}^{x}\left[B-\left|B-2\left({\mathfrak{P}}_{\mathrm{s}}^{\prime}-g0pt\right)\right|\right]\mathrm{d}x^{\prime}.$

As one can notice, surface pressure has been eliminated from the recovery
formula (17). This is because the Cauchy–Riemann system gives us two equations
(the real and complex parts) that can be combined to suppress either the
surface pressure of the surface elevation [clamond2023steady]. However, it is
still necessary to fix Bernoulli’s constant by adding a direct measure (phase
speed, wave height, etc.) or by assuming the physics at the surface (e.g.
capillary or flexural effects). For this work, we assume Stokes’s first
definition of the phase speed to be a given parameter. Thus, solving the
inverse problem reduces to finding the zeros of the set of equations (1), (12)
and (17).

###  3.4. The Levenberg–Marquardt procedure

In practice, solving the inverse problem with a direct approach is fairly
straightforward. It only requires to introduce a quadrature for numerical
integration (we choose the trapezoidal quadrature) and then insert the
different expressions in almost any root-finding algorithm (e.g. the fsolve
function in Matlab). Numerically, $\eta=\\{\eta_{i}\\}_{i}$ is a vector
defined at these quadrature points. Because this problem is notoriously
challenging to solve, we decide to use a Levenberg–Marquardt algorithm, which
has shown its efficiency before (see [C13, CC13, CH20, clamond2023steady,
CLH23, LC23]).

We quickly sketch this algorithm. Solving the inverse problem directly (for
the surface profile) can be summarised as finding the root of equation
$\mathcal{E}(\eta)=0$ for some function $\mathcal{E}$. This is equivalent as
minimising the functional
$\mathcal{J}(\eta)\vcentcolon={\textstyle{\frac{1}{2}}}\|\mathcal{E}(\eta)\|_{2}^{2}={\textstyle{\frac{1}{2}}}\mathcal{E}\bm{\cdot}\mathcal{E}$.
There are two common ways for this, both being iterative methods. Assume we
are given an initial guess $\eta^{0}$ with the overscript standing for the
iteration step. For any iteration $\ell$, we will denote
$\mathcal{E}^{\ell}\vcentcolon=\mathcal{E}(\eta^{\ell})$ and
$\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\vcentcolon=\frac{\partial\mathcal{E}}{\partial\eta}(\eta^{\ell})$
the Jacobian matrix of $\mathcal{E}$ computed at $\eta^{\ell}$.

  * -

The first method, the Gradient Descent (GD), consists in computing the
negative gradient of the function $\mathcal{J}(\eta)$ at $\eta^{\ell}$ in
order to compute $\eta^{\ell+1}$. Introducing an increment
$\delta\eta_{{}_{\text{GD}}}$, expression for the gradient follows immediately

$-\frac{\partial\mathcal{J}}{\partial\eta}(\eta^{\ell})=-\left(\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\right)^{\intercal}\mathcal{E}^{\ell}=\delta\eta_{{}_{{}_{\text{GD}}}},$

with $\mathcal{A}^{\intercal}$ denoting the transpose of a matrix
$\mathcal{A}$ and where $\partial\mathcal{E}^{\ell}/\partial\eta$ can be
approximated using finite differences. This gradient is then used in an update
formula as $\eta^{\ell+1}=\eta^{\ell}+\alpha\delta\eta_{{}_{{}_{\text{GD}}}}$,
where $\alpha$ is the length along the descent direction.

  * -

The second procedure, the Gauss–Newton (GN) method, originates from a Taylor
expansion of $\mathcal{E}$ at the first order. For some increment
$\delta\eta_{{}_{\text{GN}}}$, it yields

$\displaystyle\mathcal{J}(\eta^{\ell}+\delta\eta_{{}_{\text{GN}}})\approx$
$\displaystyle\frac{1}{2}\left(\mathcal{E}^{\ell}+\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\delta\eta_{{}_{\text{GN}}}\right)\bm{\cdot}\left(\mathcal{E}^{\ell}+\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\delta\eta_{{}_{\text{GN}}}\right)$
$\displaystyle\approx$
$\displaystyle\mathcal{J}(\eta^{\ell})+\delta\eta_{{}_{\text{GN}}}\bm{\cdot}\left(\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\right)^{\intercal}\mathcal{E}^{\ell}+\frac{1}{2}\left\|\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\delta\eta_{{}_{\text{GN}}}\right\|_{2}^{2}.$

We then minimise the function
$\delta\eta_{{}_{\text{GN}}}\mapsto\mathcal{J}(\eta^{\ell}+\delta\eta_{{}_{\text{GN}}})$
to find the optimal increment, which is equivalent to finding the root of
$\frac{\partial\mathcal{J}}{\partial\delta\eta_{{}_{\text{GN}}}}(\eta^{\ell}+\delta\eta_{{}_{\text{GN}}})=0$.
After some algebraic manipulations, it gives the linear system

(18)
$\left(\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\right)^{\intercal}\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\delta\eta_{{}_{\text{GN}}}=-\left(\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\right)^{\intercal}\mathcal{E}^{\ell}=\delta\eta_{{}_{\text{GD}}}.$

The Levenberg–Marquardt algorithm can be seen as an interpolation between the
Gradient Descent and the Gauss–Newton method. It consists in solving, instead
of (18), the Levenberg–Marquardt (LM) update as:

(19)
$\left[\left(\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}\right)^{\intercal}\frac{\partial\mathcal{E}^{\ell}}{\partial\eta}+\lambda\mathcal{I}\right]\delta\eta_{{}_{\text{LM}}}=\delta\eta_{{}_{\text{GD}}},$

where $\mathcal{I}$ is the identity matrix and $\lambda$ is a damping
parameter. A large value of $\lambda$ results in an update closer to the
Gradient Descent increment, while a small value of $\lambda$ makes the update
closer to the Gauss–Newton increment. The damping parameter is updated at each
iteration. It usually starts large and becomes slower as the iterations
converge towards the solution. This procedure is used to stabilize the
iterations when compared to the Gauss–Newton algorithm. The interested reader
can find a more thorough presentation of the Levenberg–Marquardt algorithm in
[bjorck1996numerical].

##  4\. Numerical procedure: optimisation approach

###  4.1. An optimisation formulation

The previous approach relies on the assumption that the pressure function can
be analytically continued. Without this assumption, we can formulate the
problem of recovering the surface elevation with bottom pressure measurements
within an optimisation frame. We resume the computations from (12). From
Bernoulli’s principle, we note that given $\eta$, ${p}_{\mathrm{s}}$ and $B$,
we can compute the pressure at the seabed using the formula:

(20)
${p}_{\mathrm{b}}(x)=g0pt-\frac{B}{2}-\frac{k^{2}}{32{\pi}^{2}}\left[\int_{\mathcal{C}}\cot\\!\left(k\frac{{z}_{\mathrm{s}}^{\prime}-{z}_{\mathrm{b}}}{2}\right){w}_{\mathrm{s}}^{\prime}\mathrm{d}{z}_{\mathrm{s}}-\int_{\mathcal{C}}\cot\\!\left(k\frac{{\tilde{z}}_{\mathrm{s}}^{\prime}-{z}_{\mathrm{b}}-2\mathrm{i}0pt}{2}\right){\widetilde{w}}_{\mathrm{s}}^{\prime}\mathrm{d}\widetilde{{z}_{\mathrm{s}}}^{\prime}\right]^{2}.$

This procedure is sketched in figure 1 as an input–output process.

Input $\eta$, ${p}_{\mathrm{s}}$, $B$ Black box using (20) Output
${p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)$ Figure 1. Sketch of the procedure
used for computing the bottom pressure from surface information.

Ideally, we would like to find the input leading to
${p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)={p}_{\mathrm{b}}^{\mathrm{obs}}$.
However, this strict equality may be irrelevant (due to noise in the
observation for instance) and several solutions $\eta$, ${p}_{\mathrm{s}}$ and
$B$ may exist (although not being physical). Thus, we will rather minimise
some distance between ${p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)$ and
${p}_{\mathrm{b}}^{\mathrm{obs}}$, while constraining the fields $\eta$,
${p}_{\mathrm{s}}$ and $B$ to respect the equality constraint (10) and the
mean-level condition for $\eta$. This is summarised as the following
minimisation problem

(21) $\displaystyle\min_{\eta,{p}_{\mathrm{s}},B}$
$\displaystyle\mathcal{F}(\eta,{p}_{\mathrm{s}},B)=\int_{\mathcal{P}}\left[{p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)-{p}_{\mathrm{b}}^{\mathrm{obs}}\right]^{2}$
s.t. $\displaystyle\mathcal{G}_{i}(\eta,{p}_{\mathrm{s}},B)=0,\ i=1,2,3.$

where $\mathcal{G}_{1}$, $\mathcal{G}_{2}$ and $\mathcal{G}_{3}$ correspond to
(1), (10) and (12), respectively. As mentioned, the bottom pressure
${p}_{\mathrm{b}}$ is computed from expression (20).

###  4.2. The augmented Lagrangian approach

In this section, we explain the augmented Lagrangian algorithm, that we use in
order to solve the problem (21). This is an iterative method that starts with
an initial guess $(\eta^{0},{p}_{\mathrm{s}}^{0},B^{0})$ and then, based on
the distance
$\|{p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)-{p}_{\mathrm{b}}^{\mathrm{obs}}\|_{L^{2}(\mathcal{P})}^{2}$,
will produce a new iteration $(\eta^{1},{p}_{\mathrm{s}}^{1},B^{1})$ closer to
the minimal solution. The procedure is then iterated. The main issue with (21)
are the equality constraints that should be handled with a specific procedure.
A broad range of algorithms are specifically designed for that, such as the
SQP method or the barrier methods; see [boyd2004convex] for a detailed review
on the numerical algorithms to solve constrained optimisation problems. In our
context, we opted for the augmented Lagrangian technique because of its
ability to transform constrained problems into an unconstrained one (with a
penalty term). We introduce in the following this algorithm in a generic
framework.

Given functions $\mathcal{F}$ and $\mathcal{G}$, suppose we would like to
solve the generic optimisation problem

(22) $\min_{x}\ \mathcal{F}(x)\text{ s.t. }\mathcal{G}_{i}(x)=0,\
i\in\mathcal{J}.$

where $\mathcal{J}$ is some subset of $\mathbb{N}$. The augmented Lagrangian
algorithm offers a way to handle the constraint $\mathcal{G}(x)=0$. It
consists in taking an approximate Lagrange multipliers $\lambda^{\ell}_{i}$
(in the sense of optimisation) and coefficients $\rho^{\ell}_{i}$ to solve the
optimisation problem

$x^{\ell}\in\operatorname*{arg\,min}_{x}\mathcal{F}(x)+\sum_{i\in\mathcal{J}}\lambda_{i}^{\ell}\mathcal{G}_{i}(x)+\sum_{i\in\mathcal{J}}\frac{\rho_{i}^{\ell}}{2}|\mathcal{G}_{i}(x)|^{2}.$

This procedure turns the constraint problem (22) into an unconstrained one,
with the same cost function $\mathcal{F}$ but adding the scalar product
$\lambda_{i}\mathcal{G}_{i}(x)$ and the quadratic penalization term
${\textstyle{\frac{1}{2}}}\rho_{i}^{\ell}|\mathcal{G}_{i}(x)|^{2}$. This
unconstrained minimisation problem, which will induce inner iterations, can
then be solved numerically using methods for unconstrained minimisation
problems ; we used the interior point method as implemented by the function
fmincon in Matlab (see again [boyd2004convex] for more details on this
method). Based on $x^{\ell}$ and for some tolerance $\mu^{\ell}$, the
multiplier or the coefficients are updated. For each $i\in\mathcal{J}$, it
follows that

  * -

If $|\mathcal{G}_{i}(x)|\leq\mu^{\ell}$, then
$\lambda_{i}^{\ell+1}=\lambda_{i}^{\ell}+\rho_{i}^{\ell}\mathcal{G}_{i}(x)$
and $\rho_{i}^{\ell+1}=\rho_{i}^{\ell}$ (the solution respects sufficiently
the constraints and the multiplier is thus updated).

  * -

Else, $\lambda_{i}^{\ell+1}=\lambda_{i}^{\ell}$ and
$\rho_{i}^{\ell+1}=\tau\rho_{i}^{\ell}$ for some $\tau>1$ (the constraints are
too much violated, so we increase the weight of the quadratic penalization).

An analysis of this algorithm and a proof of convergence can be found in
[bertsekas1982constrained].

(a)

(b)

(c)

(d)

Figure 2. Relative errors versus the number of Fourier modes $N$ for (a): $B$,
(b): $\eta$ and (c): ${p}_{\mathrm{s}}$. We give the running time (in hours)
in panel (d). The circled red lines correspond to the direct approach. The
crossed blue (squared magenta) lines represent the solutions from the
optimisation method with a low (high) limit on the number of function
evaluations.

##  5\. Numerical experiments

We now test numerically the two methods (direct and optimisation) on a common
example of surface profile recovery. This example was already treated in
[clamond2023steady] and appeared as one of the most difficult case considered.
We mention that every computations were made on a processor AMD EPYC 7542
@1.5Ghz with 512GB of RAM.

###  5.1. Recovery of capillary waves

When considering a fluid with surface tension $\gamma$, the analytic form of
the pressure at the surface is given by the formula

(23)
${p}_{\mathrm{s}}\vcentcolon=-\gamma\kappa=-\frac{\gamma\eta_{xx}}{\left(1+\eta_{x}^{2}\right)^{3/2}},$

where $\kappa$ is the curvature of the interface. It was shown recently
[clamond2023steady] that in this context, recovery of the surface elevation is
a challenging task due to the presence of high-order derivatives in the
nonlinear terms. Nevertheless, we demonstrate the efficiency of both
approaches by considering a gravity-capillary wave with $\gamma=1/3$,
$\sigma=-1$ and $L/d=6\pi$. The solutions of reference $B^{*}$, $\eta^{*}$ and
${p}_{\mathrm{s}}^{*}$ are computed from an algorithm adapted from [LC23] with
the total wave height being fixed at $H/0pt=0.1$. For the initial guess, we
use the solutions given by linear theory for both numerical approaches. In
addition, for the Levenberg–Marquardt algorithm, we used an initial damping
parameter set to $\lambda=0.05$. The initial parameters for the optimisation
algorithm are set to $\lambda^{0}_{i}=0$, $\rho^{0}_{i}=10$ and $\tau^{0}=2$.

First, we let vary the number of Fourier modes $N$ and observe the eventual
convergence of the solutions. In figure 2, we display with circled red lines
the relative errors on $B$, $\eta$ and ${p}_{\mathrm{s}}$, as well as the
running time, for the direct approach. We superpose the results from the
optimisation approach with a low (high) limit on the number of function
evaluations in the inner iterations using crossed blue (squared magenta)
lines. All these data are summarised in the rows of table LABEL:tab:errors
without noise. As one can notice in figure 2, the accuracy of the direct
approach continues to increase with $N$ whereas the errors from the
optimisation method slightly grow. However, restricting the number of function
evaluations on the latter allows for faster (and more accurate) results when
the number of modes is small enough. Finally, we notice that the computation
time is in favor of the direct approach by far ($\sim 100$ times faster).

(a)

(b)

(c)

Figure 3. (a): Bottom pressure of reference, computed from expression (20).
The dashed-dotted black line represents the hydrostatic value $g0pt$. (b,c):
Recovered surface pressure and surface elevation from direct approach (red
crosses) and optimisation approach (blue circles). The reference values are
displayed in dark green lines. All computations were made with $N=512$.

Eventually, we compare the recovered solutions and the reference values in
figure 3. Starting from the bottom pressure in panel 3(a), we obtain the
surface elevation and the surface pressure by direct approach (red crosses) or
by optimisation method (blue circles). Both results are in excellent agreement
with the reference value given by the dark green curves.

###  5.2. Addition of noise in the measurement

We now test the stability of the method to random perturbations of the
observation signal. For this, we add a white noise to the pressure of
reference ${p}_{\mathrm{b}}^{\mathrm{obs}}$ and try to recover solutions at
the surface and Bernoulli’s constant from this noisy data. This is of
significant importance from an engineering point of view, since most
measurements will necessarily contain noise.

Starting from the pressure field ${p}_{\mathrm{b}}^{\mathrm{obs}}$ without
noise, we create a noisy signal
${\widecheck{p}}_{\mathrm{b}}^{\mathrm{obs}}={p}_{\mathrm{b}}^{\mathrm{obs}}+\varepsilon$.
Here, the white noise $\varepsilon\sim\mathcal{N}(0,\text{STD})$ is designed
using the standard deviation
$\text{STD}=\frac{1}{3}\epsilon(\max({p}_{\mathrm{b}}^{\mathrm{obs}})-\min({p}_{\mathrm{b}}^{\mathrm{obs}}))$,
where $\epsilon$ will be some prescribed error. This let us have an error on
the signal which is at most
$\epsilon(\max({p}_{\mathrm{b}}^{\mathrm{obs}})-\min({p}_{\mathrm{b}}^{\mathrm{obs}}))$
with a 99% confidence. In our experiments, we use $\epsilon=1\%,5\%$ and
$10\%$. All results are shown in Table LABEL:tab:errors. The errors should be
understood as the relative error (with respect to the sup norm) to the
expected solutions without noise. As deduced from these results, the direct
approach fails to recover the surface even with the smallest amount of noise
added. In contrast, the optimisation process produces results in better
agreement with the solution sought. This can also be emphasised from figure
5.2, where we used $N=128$ modes and $\epsilon=5\%$. We notice that all errors
are of the same magnitude as the noise added to the original signal, exception
made of ${p}_{\mathrm{s}}$. However, the computation time needed to solve the
inverse problem also grows significantly although it is not proportional to
the rate of noise added.

% noise | N |  Direct |  Optim.
---|---|---|---
[ head to column names, filter expr = test  (a) |  |  |
|  |  |

##  6\. Discussion

We now discuss the results presented above and draw some conclusions from
them. We will mainly discuss two aspects: the effect of noise addition on the
two methods, and the comparison of the direct and optimisation methods from a
technical perspective. Both topics will be the occasion to hint some possible
extensions and future research directions.

###  6.1. Comments on the noise addition

Implementing noise in our numerical experiments provided some evidence that
analytical continuation of the bottom pressure is a prerequisite for the
direct method to work. Even in the case with only 1% of error, results in
table LABEL:tab:errors show that it renders this method dramatically
ineffective. In contrast, the optimisation procedure shows a good agreement
even in the presence of noise. For instance, the solution fields $\eta$ and
$B$ have a relative error of the same magnitude as the noise added, which
demonstrates the robustness of this method. Actually, this observation seems
natural. The way we formulate the inverse problem using an optimisation frame
can be seen as a filtering for the bottom pressure, since minimizing the
$L_{2}$ distance between ${p}_{\mathrm{b}}(\eta,{p}_{\mathrm{s}},B)$ and
${\widecheck{p}}_{\mathrm{b}}^{\mathrm{obs}}$ corresponds to fitting
${p}_{\mathrm{b}}$ to some unbiased estimator of
${\widecheck{p}}_{\mathrm{b}}^{\mathrm{obs}}$. However, the surface pressure
${p}_{\mathrm{s}}$ still shows a non-negligible deviation from the sought
solution. As mentioned in the introduction, this can be attributed to the ill-
posedness of the problem, namely that any disturbances in the measurements
propagate exponentially to the surface. This is further explained by the
nature of the dynamic condition at the surface and the physical assumptions on
${p}_{\mathrm{s}}$. Indeed, capillary (or flexural) waves involve high-order
derivatives of the surface profile that are numerically problematic to
approximate when the solution is not smooth enough (as it is the case with
unfiltered noise). Numerical errors can thus accumulate since there are no
constraint to enforce any regularity on the surface pressure. We see several
ways to accommodate both methods to the case where noise is added. Regarding
the direct approach, we should focus on filtering the noise in the bottom
measurements before solving the analytical expressions. This method would
require a precise estimation technique to guarantee the analyticity of the
signal, since our experiments show that even the slightest addition of noise
breaks the whole recovery process. As for the optimisation method, the main
issue in the presence of noise solely concerns the surface pressure
${p}_{\mathrm{s}}$. This can be solved by including more information to the
problem, e.g. some additional measurements or some knowledge of the physical
mechanisms (the two being not exclusive). In practice, this added information
is simply incorporated in the optimisation problem by adding an extra
constraint. However, the latter could make the problem unsolvable or at least
numerically harder to solve. We illustrate this idea with a numerical
experiment on the same configuration than before where we added the constraint
$\|{p}_{\mathrm{s}}-{p}_{\mathrm{s}}^{\textrm{th}}(\eta)\|^{2}_{L^{2}(\mathcal{P})}=0$
in the optimisation problem (21) (${p}_{\mathrm{s}}^{\textrm{th}}$ being
computed from the analytic expression (23)). The results for $N=128$ and
$\epsilon=5\%$ are depicted in figure 5.2 with a dashed red curve when adding
this new constraint and in dotted red line without it. As clear from these
panels, knowing the physics at the surface removes the oscillations due to the
noise and allows to recover the solutions of reference with an excellent
agreement. However, this configuration becomes extremely challenging to solve
when considering this additional constraint, as clear from the computation
time taken to solve this problem (which was $t\sim 5.1243$ hours). Future
research should focus on the nature of information one can add to stabilise
the solutions and lower the complexity.

###  6.2. Comparison of the direct and optimisation approaches

The principal objective of this work was to solve a nonlinear ill-posed inverse problem using two distinct approaches and compare these methods from a numerical perspective. At first glance, we notice a clear discrepancy in the computation time of both procedures. The direct approach is indeed much faster than the optimisation method and we would like to understand the reasons behind this. Let us return to the intrinsic nature of our direct formulation. By considering the analytic extension of the bottom pressure, it allows for the derivation of some implicit formula for the recovery of the surface profile. Solving these expressions with a root-finding algorithm represents an efficient way of solving this inverse problem with a system of reduced size. Therefore, the direct approach naturally exhibits one of the fastest way to compute these variables with an excellent accuracy. On the other hand, the optimisation procedure (21) requires to solve a series of unconstrained problem iteratively and adjust the parameters at each iteration to respect the constraints. This method acts on a system with a larger size and is strongly dependent on the number of constraints and optimisation variables. Indeed, we observe in figure 2 (with data summarised in table LABEL:tab:errors with $0\%$ error) a tendency of the errors on the recovered variables to grow as the size of the system increases, while the direct method takes advantage of this increased number of variables. This is most probably due to the complexity of the formulation (i.e. respecting all the constraints simultaneously) becoming too arduous. Acceleration of the optimisation method should be an important topic of research in order to make this approach really efficient and could be made using different techniques, such as an increased parallelisation of the computation or via techniques of distributed optimisation ; see [bertsekas2015parallel, boyd2011distributed] for more details on this topic. Moreover, both methods are interesting in the sense that they can be generalised fairly easily to more complex situations. Adding a linear shear current (as done in [CLH23]) and considering overhanging profiles (as done in [LC23]) or flexural effects at the surface is straightforward to implement in both approaches, although it would increase the complexity of the expressions we used. Extending this work to the case of unsteady motions is a challenging task that is still lacking from the literature. It would require in practice to solve the inverse problem iteratively in time (repeating the current technique at each time step) while avoiding numerical instabilities (e.g. aliasing) and following the same branch of solution. In principle, both methods are applicable when considering this configuration although the running time of the optimisation problem represents a clear impediment to this formulation. Finally, extension towards the three-dimensional case is most certainly the principal reason why we considered the optimisation approach in the first place. This formulation is indeed well-designed to handle three-dimensional problems ; see [challis2009level, gaymann2019fluid, villanueva2017cutfem] for examples of optimisation problems involving 3D equations. However, it would require to solve the equations of motion directly instead of considering the analytical expressions we used in this work that were obtained from a boundary integral method. Since these expressions are expressed in the complex plane, their generalisation in 3D is not trivial and represents a current topic of research. For this reason, the direct approach is less inclined to be extended to this configuration whereas the optimisation formulation could present a solution to this long-standing problem. |  |  |
Funding. Joris Labarbe has been supported by the French government, through the $\mbox{UCA}^{\mbox{\tiny JEDI}}$ Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01. |  |  |

(b)

(c)

(a) (a): Pressure at the seabed. Ground truth value (computed from (20)) is
depicted in blue line. The signal with $5\%$ added noise is represented in
dotted red line. (b,c): Recovered surface pressure and surface elevation
without (with) knowledge of the physics at the surface are represented in
dotted (dashed) red lines. The reference values are displayed with blue lines.
All computations were made with $N=128$.