# Identification of bottom deformations of the ocean from surface measurements

Laurent Bourgeois, Jean-François Mercier, Raphaël Terrine

► To cite this version:

Laurent Bourgeois, Jean-François Mercier, Raphaël Terrine. Identification of bottom deformations of the ocean from surface measurements. Inverse Problems and Imaging, 2025, 19 (6), pp.1075-1113. 10.3934/ipi.2025008. hal-04514543v2

HAL Id: hal-04514543

https://inria.hal.science/hal-04514543v2

Submitted on 28 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Identification of bottom deformations of the ocean
from surface measurements

January 28, 2025

LAURENT BOURGEOIS, JEAN-FRANÇOIS MERCIER AND RAPHAËL TERRINE

Laboratoire Poems, CNRS, INRIA, Ensta Paris, Institut Polytechnique de Paris,
828 Boulevard des Maréchaux, 91762 Palaiseau, France

## Abstract

In this paper we consider a general scheme to solve two different inverse problems related to oceanography, that is retrieving either a tsunami or the shape of the seabed from the measurement of the free surface perturbation. We consider two dimensional geometries and linear potential models in the frequency regime. Such general scheme consists in firstly recovering the potential in the whole domain and secondly compute the seeked parameter at the bottom of the ocean, which in the two inverse problems is a function involved in a more or less complicated boundary condition. The first step amounts to solve an ill-posed Cauchy problem for the Laplace or Helmholtz equation, which we regularize by using a mixed formulation of the Tikhonov regularization and the Morozov principle to compute the regularization parameter. The computation of such Tikhonov-Morozov solution is based on an iterative method consisting in solving a sequence of weak formulations which are discretized with the help of a simple Lagrange type Finite Element Method. In the particular case of the acoustic model, we need to solve a Laplace-type equation associated with the noisy Neumann boundary data and compute the noise amplitude of its solution. A probabilistic method is proposed to obtain such amplitude of noise. Some numerical experiments show the feasibility of our strategy.

## 1 Introduction

This paper is devoted to two different inverse problems which arise in oceanography. The first one is the identification of a tsunami from measurements of the free surface deformation, such tsunami being characterized by a brutal displacement of the bottom surface of the ocean. The second one is the bathymetry problem, which consists in recovering the underwater depth by
using the same measurements as in the previous problem. In these two problems, the goal is to retrieve some sea bottom parameters from surface data, but while the loading is passive in the first problem (the tsunami is a natural phenomenon), it is active in the second one (a source artificially generated is required).

Concerning the literature on tsunamis, there is a large amount of contributions related to the forward modeling of tsunamis, for example [33, 34, 36]. In particular, a comprehensive study of the behaviour of a tsunami at all steps of its “life” and of the best models to use at each of these steps is exposed in [18]. However, the identification of tsunamis from surface measurements is far less addressed in the literature. In this direction, let us mention [26], in which a shallow-water model is handled and a least-square method is proposed, and [24], which is based on a machine-learning algorithm. For the simpler linear water wave problem in the presence of a flat bottom, in [14] a numerical approach based on integral transforms with respect to the horizontal directions and with respect to time is successfully used to invert experimental data. This paves the way to real applications of inverse problems in oceanography. The literature on bathymetry seen as an inverse problem is also quite poor. Let us mention [35], in which a numerical reconstruction procedure is conducted, in the presence of a nonlinear water wave model and by using periodic boundary conditions in the horizontal directions. We also mention [22, 32], which address the theoretical question of uniqueness/stability for the inverse problem in the context of the nonlinear water wave problem. Note that the well-posedness of the corresponding forward problem is proved in [29].

In this contribution we address the two previous inverse problems by restricting ourselves to some two-dimensional linear models in the frequency domain. The physics of the problem is then fully governed by the scalar velocity potential. Note that such linear models are only valid during the first step of the tsunami generation and would not be appropriate during the second propagation step, as emphasized in [18]. The two inverse problems, which are strongly ill-posed, are attacked by using the same strategy, which consists in two steps. In the first step, surface data are used to compute an approximate velocity potential in the whole volumic domain by transforming the initial ill-posed problem into a well-posed regularized one. In the second step we use the partial derivatives of such potential on the bottom of the domain to compute either the time derivative of the vertical displacement which characterizes the tsunami (first inverse problem) or the geometric perturbation of the bottom shape with respect to the horizontal flat reference (second inverse problem). Note that in the bathymetry problem, we are led to make a strong assumption on the smoothness of the seabed to ensure uniqueness for the inverse problem (see Theorem 4.2). The novelty of our contribution is twofold. Firstly in these two steps we use weak formulations
and finite element methods, which are well adapted to non-flat sea bottoms. We observe that when the bottom is flat, the forward tsunami problem is often solved in the literature by first using a Fourier transform in the horizontal directions, which is forbidden in the presence of a variable sea bottom. Secondly we wish to study the influence of the choice of the model on the quality of the reconstruction at different frequencies. In this view we start from a model which takes account of both acoustic and gravity waves, the former being dominating at high frequencies, the latter being dominating at low frequencies.

The main ingredient of the first step described above is the notion of mixed formulation of the Tikhonov regularization, which is designed to regularize ill-posed Cauchy problems and can be immediately discretized with the help of a classical Lagrange Finite Element Method. This technique was first introduced in [5] in the case of the Cauchy problem for the Laplace equation, then generalized to an abstract context in [12]. The Tikhonov formulation involves a small regularization parameter that we choose following the Morozov principle (see for example [28]) and compute numerically by using the Fenchel-Rockafellar duality in optimization (see the general theory in [19]). This technique was first proposed in the context of the Morozov principle in [6, 10] again in the case of the Cauchy problem for the Laplace equation, then generalized to an abstract setting in [11]. The duality approach leads to the problem of minimizing a convex but non-quadratic functional, which in practice turns out to be delicate. This is why, following an idea of [16], we transform such problem into a sequence of quadratic minimization problems, the solutions of which converge to the solution of the original problem. Such quadratization technique has already been used in [11].

Our paper is organized as follows. In section 2 we present the linear model of oceanography that we will consider, that is the so-called complete model, and we address the corresponding forward problems, both for the tsunami problem and the bathymetry problem. Two classical limit models deriving from the complete model, that is the gravity model and the acoustic model, are presented in section 3. Our two inverse problems, that is the tsunami identification and the bathymetry problems, are introduced in section 4. Section 5 is devoted to the abstract framework of the mixed formulation of the Tikhonov regularization method, such framework being the core of our strategy to solve these two inverse problems, whatever the chosen model is. This abstract framework is then applied to the tsunami identification problem and to the bathymetry problem in sections 6 and 7, respectively. To solve the bathymetry problem, assuming the geometric perturbation is small, we consider a second order impedance boundary condition on a flat bottom instead of the true perturbed bottom, which is an approximation based on an asymptotic analysis. Note that we assume a lower bound of
the surface derivative of the approximate velocity potential at the bottom (see Assumption 7.2). In the case of the acoustic model, we need to convert the amplitude of noise that contaminates a Neumann boundary data to the amplitude of noise of a volumic lifting solution of such Neumann boundary problem. To proceed, two approaches, a deterministic and a probabilistic one, are exposed in section 8. Such section is complemented by an appendix, in which we provide more theoretical results related to the probabilistic approach in a slightly simpler case from the point of view of the geometry. Lastly, numerical experiments are presented in section 9.

## 2 The selected model and the forward problems

In this article, the ocean at rest is defined as a two-dimensional domain $\Omega_\theta$, given by

$$ \Omega_\theta = \{(x, z) \in \mathbb{R}^2, -H + \theta(x) < z < 0\}, $$

where $\theta$ is a continuous and compactly supported function which can be considered as a perturbation of the flat bottom of equation $z = -H$, while the free surface at rest is characterized by $z = 0$. The domain $\Omega_\theta$ is hence delimited by the two curves

$$ \Gamma_0 = \{(x, z) \in \mathbb{R}^2, z = 0\}, \quad \Gamma_{-H+\theta} = \{(x, z) \in \mathbb{R}^2, z = -H + \theta(x)\}. $$

### 2.1 The complete model

We first consider the tsunami problem. Let us denote $\zeta(x, t)$ the upward normal displacement of the sea bottom, which is the only source of the problem, and $s(x, t)$ the induced perturbation of the free surface. The velocity potential $\varphi(x, z, t)$, the functions $s$ and $\zeta$ are related to each other by the following equations and boundary conditions:

$$ \left\{ \begin{array}{l} (1/c^2)\partial_t^2\varphi - \Delta\varphi = 0 \quad \text{in } \Omega_\theta \times (0, +\infty), \\ \partial_t\varphi + gs = 0 \quad \text{on } \Gamma_0 \times (0, +\infty), \\ \partial_z\varphi - \partial_t s = 0 \quad \text{on } \Gamma_0 \times (0, +\infty), \\ \partial_\nu\varphi = -\partial_t\zeta \quad \text{on } \Gamma_{-H+\theta} \times (0, +\infty) \end{array} \right. \quad (1) $$

where $\nu$ is the unit outward normal to $\Omega_\theta$, $g$ is the gravitational acceleration and $c$ is the celerity of the acoustic waves. The system (1) has to be complemented by initial conditions. Note that the two boundary conditions on $\Gamma_0$ can be merged into the single equation:

$$ \partial_t^2 \varphi + g \partial_z \varphi = 0. \quad (2) $$

The linear model (1) is rigorously derived from the Euler equations in [17] and corresponds to the so-called “barotropic case” in [17], with the additional assumption that the density is constant.
In the frequency domain, which amounts to assume that $\varphi(x, z, t) = u(x, z)e^{-i\omega t}$, $s(x, t) = \eta(x)e^{-i\omega t}$ and $\zeta(x, t) = \chi(x)e^{-i\omega t}$ for some fixed frequency $\omega > 0$, the forward scattering problem of the tsunami reads: for $\chi \in L^2(\Gamma_0)$, where $\chi$ is supposed to be compactly supported, find $u \in H_{\text{loc}}^1(\Omega_\theta)$ such that

$$ \left\{ \begin{array}{lll} \Delta u + (\omega^2/c^2)u = & 0 & \text{in } \Omega_\theta, \\ \partial_z u - (\omega^2/g)u = & 0 & \text{on } \Gamma_0, \\ \partial_\nu u = & i\omega\chi & \text{on } \Gamma_{-H+\theta}, \\ u & \text{is} & \text{outgoing.} \end{array} \right. \qquad (3) $$

Here,

$$ H_{\text{loc}}^1(\Omega_\theta) = \{u \in \mathcal{D}'(\Omega_\theta), \psi(x)u \in H^1(\Omega_\theta), \forall \psi \in C_0^\infty(\mathbb{R})\}. $$

The last line of (3) corresponds to a radiation condition, which is necessary to guarantee uniqueness in a scattering problem, and which will be clarified hereafter. Note that once the potential $u$ is known, the perturbation of the free surface is given by either of the two formulas

$$ \eta = \frac{i\omega}{g}u, \quad \eta = -\frac{1}{i\omega}\partial_z u. $$

The system (1) is interesting in the context of tsunamis because it relies on a model which involves both acoustic and gravity waves, in other words both constants $g$ and $c$ appear in (1) (see for example [33, 34, 36]). Such model will be called the complete model in the following.

Now let us introduce the forward bathymetry problem: for some source point $N \in \Omega_\theta$, find $u \in H_{\text{loc}}^1(\Omega_\theta)$ such that

$$ \left\{ \begin{array}{lll} -\Delta u - (\omega^2/c^2)u = & \delta_N & \text{in } \Omega_\theta, \\ \partial_z u - (\omega^2/g)u = & 0 & \text{on } \Gamma_0, \\ \partial_\nu u = & 0 & \text{on } \Gamma_{-H+\theta}, \\ u & \text{is} & \text{outgoing.} \end{array} \right. \qquad (4) $$

Note that in problem (4), the source term appears in the right-hand side of the first equation, and no more in the boundary condition on the seabed $\Gamma_{-H+\theta}$. Let us prove that the forward problems for the complete model, both in the case of the tsunami generation (3) and the bathymetry (4), are well-posed. In this view, we have to specify what we mean by “outgoing” solutions in the above forward problems, which leads us to introduce the modes.

## 2.2 Determination of the modes

Like in the pioneering work [30], we introduce the modes, that is the solutions in $\Omega_0 = \mathbb{R} \times (-H, 0)$ in the form $u(x, z) = f(z)e^{\lambda x}$ to the following problem
in the absence of the source term $\chi$ or $\delta_N$:

$$ \left\{ \begin{array}{l} \Delta u + (\omega^2/c^2)u = 0 \quad \text{in } \Omega_0, \\ \partial_z u - (\omega^2/g)u = 0 \quad \text{on } \Gamma_0, \\ \partial_\nu u = 0 \quad \text{on } \Gamma_{-H}. \end{array} \right. $$

Finding such $u$ amounts to find $(\lambda, f) \in \mathbb{C} \times H^1(-H, 0)$, $f \neq 0$, such that

$$ \left\{ \begin{array}{l} f''(z) + (\lambda^2 + \omega^2/c^2)f(z) = 0 \quad \text{in } (-H, 0), \\ f'(0) - (\omega^2/g)f(0) = 0, \\ f'(-H) = 0. \end{array} \right. $$

Let us denote $\nu_0$ the unique positive solution to the equation

$$ \nu \tanh(\nu H) = \omega^2/g, \qquad (5) $$

while the sequence $(\nu_n)_{n \ge 1}$ is formed by the positive and increasing solutions to the equation

$$ \nu \tan(\nu H) = -\omega^2/g. \qquad (6) $$

Let us in addition assume the following.

**Assumption 2.1.** The frequency $\omega$ is such that there are no $n \in \mathbb{N}$, $n \ge 1$, satisfying

$$ \nu_n = \omega/c. $$

Such assumption implies that there exists some $P \in \mathbb{N}$, such that

$$ \left\{ \begin{array}{ll} \nu_n < \omega/c & \text{if } 1 \le n \le P, \\ \nu_n > \omega/c & \text{if } n > P. \end{array} \right. $$

The solutions $(\lambda, f)$ to the above problem are given by:

$$ \left\{ \begin{array}{ll} \pm i \sqrt{\nu_0^2 + \omega^2/c^2}, & f_0(z) := A_0 \cosh(\nu_0(z + H)), \\ \pm i \sqrt{\omega^2/c^2 - \nu_n^2}, & f_n(z) := A_n \cos(\nu_n(z + H)), \quad 1 \le n \le P, \\ \pm \sqrt{\nu_n^2 - \omega^2/c^2}, & f_n(z) := A_n \cos(\nu_n(z + H)), \quad n > P, \end{array} \right. \qquad (7) $$

where the $A_n$ are arbitrary non zero constants. The modes are then given by the functions

$$ u_n^{\pm}(x, z) = f_n(z)e^{\pm\beta_n x}, \qquad (8) $$

where from now on the constants $A_n$ of (7) are chosen such that the $(f_n)$ form an orthonormal complete basis of $L^2(-H, 0)$ (which is possible from a classical result of spectral theory for self-adjoint operators) and

$$ \left\{ \begin{array}{l} \beta_0 := i \sqrt{\nu_0^2 + \omega^2/c^2}, \\ \beta_n := i \sqrt{\omega^2/c^2 - \nu_n^2}, \quad 1 \le n \le P, \\ \beta_n := -\sqrt{\nu_n^2 - \omega^2/c^2}, \quad n > P. \end{array} \right. \qquad (9) $$
The modes $u_n^+$ are either propagating (for $0 \le n \le P$) or evanescent (for $n > P$) in the right direction of $\Omega_0$, so that there are outgoing to the right. The modes $u_n^-$ are either propagating (for $0 \le n \le P$) or evanescent (for $n > P$) in the left direction of $\Omega_0$, so that there are outgoing to the left.

We are now in a position to define an outgoing function $u$ in $\Omega_\theta$ as a function for which there exists $R > 0$ such that

$$u(x, z) = \sum_{n \in \mathbb{N}} c_n^{\pm} u_n^{\pm}(x, z), \quad \pm x > R, \qquad (10)$$

where the $c_n^{\pm}$, $n \in \mathbb{N}$, are complex numbers. More precisely, $R$ is any sufficiently large real such that the support of $\theta$ and $\chi$ are contained in $(-R, R)$. Let us consider the problem (3) in the simplified case of the flat seabed, that is $\theta = 0$. It will be always the case for the tsunami problem.

## 2.3 Well-posedness of the tsunami problem

In order to prove the well-posedness of (3), we introduce the bounded domain $\Omega_0^R := \{(x, z) \in \Omega_0, -R < x < R\}$ delimited by the transverse sections $\Sigma_{\pm R} = \{\pm R\} \times (-H, 0)$, assuming that the support of $\chi$ is contained in $(-R, R)$. Classically, the problem (3) in the unbounded domain $\Omega_0$ is equivalent to the following problem set in the bounded domain $\Omega_0^R$: find $u \in H^1(\Omega_0^R)$ such that

$$\left\{
\begin{array}{l}
\Delta u + (\omega^2/c^2)u = 0 \quad \text{in } \Omega_0^R, \\
\partial_z u - (\omega^2/g)u = 0 \quad \text{on } \Gamma_0^R, \\
\partial_z u = -i\omega\chi \quad \text{on } \Gamma_{-H}^R, \\
\pm\partial_x u - T_{\pm}^c u = 0 \quad \text{on } \Sigma_{\pm R},
\end{array}
\right.
\qquad (11)$$

where $\Gamma_0^R$ and $\Gamma_{-H}^R$ are the subparts of $\Gamma_0$ and $\Gamma_{-H}$ delimited by the transverse sections $\Sigma_{-R}$ and $\Sigma_R$, respectively, while the Dirichlet-To-Neumann operators $T_{\pm}^c$ are defined by

$$\left\{
\begin{array}{l}
T_{\pm}^{c}: H^{1/2}(\Sigma_{\pm R}) \rightarrow \tilde{H}^{-1/2}(\Sigma_{\pm R}) \\
\varphi \mapsto \sum_{n \in \mathbb{N}} \beta_n(\varphi, f_n)_{L^2(-H,0)} f_n.
\end{array}
\right.
\qquad (12)$$

Here, $\tilde{H}^{-1/2}(\Sigma_{\pm R})$ is the dual space of $H^{1/2}(\Sigma_{\pm R})$, the functions $f_n$ are given by (7) and the complex numbers $\beta_n$ are given by (9). It is easy to check that in problem (11), the artificial boundary conditions on the transverse sections $\Sigma_{\pm R}$ are equivalent to the radiation condition (10). For a proof in a slightly more difficult case, the reader may refer to [4].

An equivalent weak formulation to the strong problem (11) is: find $u \in H^1(\Omega_0^R)$ such that for all $v \in H^1(\Omega_0^R)$,

$$a(u, v) = l(v), \qquad (13)$$
where the sesquilinear form $a$ and the antilinear form $\ell$ are given by

$$ 
\begin{aligned}
a(u, v) &= \int_{\Omega_0^R} (\nabla u \cdot \nabla \bar{v} - (\omega^2/c^2)u \bar{v}) \, dxdz \\
&\quad - (\omega^2/g) \int_{\Gamma_0^R} u \bar{v} \, dx - \langle T_+^c u, \bar{v} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\
&\quad - \langle T_-^c u, \bar{v} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})},
\end{aligned} 
\quad (14) $$

where $\langle \cdot, \cdot \rangle_{H^{-1/2}(\Sigma_{\pm R}), \tilde{H}^{1/2}(\Sigma_{\pm R})}$ denotes the duality bracket between $H^{-1/2}(\Sigma_{\pm R})$ and $\tilde{H}^{1/2}(\Sigma_{\pm R})$ and

$$ \ell(v) = i\omega \int_{\Gamma_{-H}^R} \chi \bar{v} \, dx. \quad (15) $$

We have the following well-posedness result.

**Theorem 2.2.** The weak formulation (13) has a unique solution.

*Proof.* We introduce the decomposition $a = b + c$, where the sesquilinear forms $b$ and $c$ are given by

$$ 
\left\{ 
\begin{aligned}
b(u, v) &= \int_{\Omega_0^R} (\nabla u \cdot \nabla \bar{v} + u \bar{v}) \, dxdz - \langle T_+^c u, \bar{v} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\
&\quad - \langle T_-^c u, \bar{v} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})}, \\
c(u, v) &= -(1 + (\omega^2/c^2)) \int_{\Omega_0^R} u \bar{v} \, dxdz - (\omega^2/g) \int_{\Gamma_0^R} u \bar{v} \, dx,
\end{aligned}
\right. 
$$

and define the operators $B, C: H^1(\Omega_0^R) \to H^1(\Omega_0^R)$ such that

$$ (Bu, v)_{H^1(\Omega_0^R)} = b(u, v), \quad (Cu, v)_{H^1(\Omega_0^R)} = c(u, v), \quad \forall u, v \in H^1(\Omega_0^R). $$

In view of (9), we have that

$$ \mathcal{R}e\{b(u, u)\} = \|u\|_{H^1(\Omega_0^R)}^2 + \sum_{\pm, n>P} \sqrt{\nu_n^2 - \omega^2/c^2} |(u|_{\Sigma_{\pm R}}, f_n)_{L^2(-H,0)}|^2 \ge \|u\|_{H^1(\Omega_0^R)}^2, $$

where $\mathcal{R}e\{\cdot\}$ denotes the real part. Then $B$ is an isomorphism while $C$ is a compact operator, the Fredholm alternative allows us to conclude that in problem (13), existence is equivalent to uniqueness. In order to prove uniqueness (having in mind that $\theta = 0$), assume that $u$ satisfies the problem (11) for $\chi = 0$. From the three first equations of (11), we obtain that $u$ is an infinite linear combination of the modes $u_n^\pm$ given by (8). Since in addition the solution $u$ is outgoing, that is both leftgoing and rightgoing as specified by (10), such linear combination is necessarily equal to 0. This complete the proof. $\square$
## 2.4 Well-posedness of the bathymetry problem

Now let us prove well-posedness of the problem (4). In this view, we introduce the fundamental solution in the waveguide $\Omega_0$, that is for some source point $N(x', z') \in \Omega_0$, the solution $\mathcal{G}^c(\cdot, N) \in L_{\text{loc}}^2(\Omega_0)$ to the problem

$$ egin{cases} -\Delta \mathcal{G}^c(\cdot, N) - (\omega^2/c^2)\mathcal{G}^c(\cdot, N) = \delta_N & \text{in } \Omega_0, \\ \partial_z \mathcal{G}^c(\cdot, N) - (\omega^2/g)\mathcal{G}^c(\cdot, N) = 0 & \text{on } \Gamma_0, \\ \partial_z \mathcal{G}^c(\cdot, N) = 0 & \text{on } \Gamma_{-H}, \\ \mathcal{G}^c(\cdot, N) & \text{is outgoing.} \end{cases} $$

A short computation shows that an explicit expression of $\mathcal{G}^c$ is given, for any $M(x, z)$, by

$$ \mathcal{G}^c(M, N) = - \sum_{n \in \mathbb{N}} \frac{1}{2\beta_n} e^{\beta_n |x-x'|} f_n(z) f_n(z'), \quad (16) $$

where the $f_n$ and the $\beta_n$ are given by (7) and (9), respectively. In what follows we will assume that the support of $\theta$ lies between $\Sigma_{-R}$ and $\Sigma_R$ and that $N$ is located such that $x' < -R$. By defining the incident field $u^i$ in $\Omega_\theta$ by

$$ u^i = \begin{cases} \mathcal{G}^c(\cdot, N) & \text{if } x < -R, \\ 0 & \text{if } x > -R, \end{cases} \quad (17) $$

and by writing the problem satisfied by the corresponding outgoing scattered field $u^s = u - u^i$, we obtain that the solution $u$ to the problem (4) also satisfies in $\Omega_\theta^R$:

$$ egin{cases} \Delta u + (\omega^2/c^2)u = 0 & \text{in } \Omega_\theta^R, \\ \partial_z u - (\omega^2/g)u = 0 & \text{on } \Gamma_0^R, \\ \partial_\nu u = 0 & \text{on } \Gamma_{-H+\theta}^R, \\ -\partial_x u = T_-^c u - 2\partial_x \mathcal{G}^c(\cdot, N) & \text{on } \Sigma_{-R}, \\ \partial_x u = T_+^c u & \text{on } \Sigma_R. \end{cases} \quad (18) $$

**Remark 1.** We note that the fundamental solution $\mathcal{G}^c(\cdot, N)$ does not satisfy a Neumann boundary condition on the domain $\Omega_\theta$ (only on $\Omega_0$). This explains why we have introduced a piecewise incident field $u^i$ in (17). The fact that the incident field $u^i$ and the scattered field $u^s$ have a discontinuity across the section $x = -R$ may seem surprising but is mathematically correct, as long as the total field is continuous across such section.

An equivalent weak formulation to the strong problem (18) is: find $u \in H^1(\Omega_\theta^R)$ such that for all $v \in H^1(\Omega_\theta^R)$,

$$ a_\theta(u, v) = l_\theta(v), \quad (19) $$
where the sesquilinear form $a$ and the antilinear form $\ell$ are given by

$$ a_{\theta}(u, v) = \int_{\Omega_{\theta}^{R}} (\nabla u \cdot \nabla \bar{v} - (\omega^2/c^2)u \bar{v}) dxdz - (\omega^2/g) \int_{\Gamma_{0}^{R}} u \bar{v} dx \\ - \langle T_{+}^{c} u, \bar{v} \rangle_{H^{-1/2}(\Sigma_{R}), \tilde{H}^{1/2}(\Sigma_{R})} - \langle T_{-}^{c} u, \bar{v} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \quad (20) $$

and

$$ \ell_{\theta}(v) = -2 \int_{\Sigma_{-R}} \partial_{x} \mathcal{G}^{c}(\cdot, N) \bar{v} \, dx. \quad (21) $$

We have the following well-posedness result in the Fredholm sense.

**Theorem 2.3.** Assume that the weak formulation (19) has at most one solution. Then such weak formulation is well-posed.

The proof of Theorem 2.3 relies on a decomposition of the sesquilinear form $a_{\theta}$ of the same type as in the proof of Theorem 2.2. The difference is that uniqueness is not guarranteed. Indeed, due to the fact that the waveguide $\Omega_{\theta}$ is not straight, some trapped solutions may exist.

# 3 Two limit cases

## 3.1 Derivation of the two limits

Two well-known limit cases can be derived from the complete model. In the case of the tsunami problem for example, starting from problem (3), the gravity model (also often called the water wave model) corresponds to the case when $c \to +\infty$, that is the problem (3) becomes: find $u \in H_{\text{loc}}^{1}(\Omega_{\theta})$ such that

$$ \left\{ \begin{array}{l} \Delta u = 0 \quad \text{in } \Omega_{\theta}, \\ \partial_{z} u - (\omega^{2}/g)u = 0 \quad \text{on } \Gamma_{0}, \\ \partial_{\nu} u = i\omega\chi \quad \text{on } \Gamma_{-H+\theta}, \\ u \text{ is outgoing.} \end{array} \right. \quad (22) $$

The acoustic model corresponds to the case when $g \to 0$, that is the problem (3) becomes: find $u \in H_{\text{loc}}^{1}(\Omega_{\theta})$ such that

$$ \left\{ \begin{array}{l} \Delta u + (\omega^2/c^2)u = 0 \quad \text{in } \Omega_\theta, \\ u = 0 \quad \text{on } \Gamma_0, \\ \partial_\nu u = i\omega\chi \quad \text{on } \Gamma_{-H+\theta}, \\ u \text{ is outgoing.} \end{array} \right. \quad (23) $$

One can wonder which model we have to choose. To decide, we can define two typical frequencies $\omega_g := \sqrt{g/H}$ and $\omega_a := c/H$, such that for typical realistic values of $g, c, H$ we have $\omega_a \gg \omega_g$. The gravity model (22) is well adapted for $\omega \sim \omega_g$ while the acoustic model (23) is well adapted for $\omega \sim \omega_a$.
## 3.2 Analysis of the gravity case model

The results obtained for the complete model in the previous section are valid in the gravity case, except that when $c \to +\infty$ the modes have a simpler form. They are given by (8), where in the definition of the $f_n$ given by (7) and the definition of the $\beta_n$ given by (9) we have $P = 0$, which yields

$$ \begin{cases} f_0(z) := A_0 \cosh(\nu_0(z + H)), \\ f_n(z) := A_n \cos(\nu_n(z + H)), \quad n \ge 1, \end{cases} \qquad (24) $$

and

$$ \begin{cases} \beta_0 := i\nu_0, \\ \beta_n := -\nu_n, \quad n \ge 1. \end{cases} \qquad (25) $$

In the gravity case, to avoid ambiguity the Dirichlet-To-Neumann operators given by (12) will be denoted $T_{\pm}^g$ in the sequel.

## 3.3 Analysis of the acoustic model

Let us now address the forward acoustic problems, both in the case of the tsunami generation and the bathymetry. The analysis is slightly different from the complete model. Following the same lines as in section 2, we need again to compute the modes, that is the solutions in $\Omega_0 = \mathbb{R} \times (-H, 0)$ in the form $u(x, z) = g(z)e^{\mu x}$ to the problem

$$ \left\{ \begin{array}{ll} \Delta u + (\omega^2/c^2)u = 0 & \text{in } \Omega_0, \\ u = 0 & \text{on } \Gamma_0, \\ \partial_z u = 0 & \text{on } \Gamma_{-H}. \end{array} \right. \qquad (26) $$

The derivation of these modes from the modes of the complete model is not obvious, this is why we compute them starting from (26). What follows relies on the following assumption.

**Assumption 3.1.** The frequency $\omega$ is such that there are no $n \in \mathbb{N}$ satisfying

$$ \mu_n = \frac{\omega}{c}, \quad \mu_n := \left(n + \frac{1}{2}\right) \frac{\pi}{H}. $$

Such assumption implies that there exists some $P \in \mathbb{N}$, $P > 0$, such that

$$ \begin{cases} \mu_n < \omega/c & \text{if } 0 \le n < P, \\ \mu_n > \omega/c & \text{if } n \ge P. \end{cases} $$

If we denote

$$ g_n(z) = \sqrt{2/H} \cos(\mu_n(z + H)), \quad n \in \mathbb{N}, \qquad (27) $$

which implies that $\|g_n\|_{L^2(-H,0)} = 1$, and

$$ \left\{ \begin{array}{ll} \gamma_n = i\sqrt{(\omega^2/c^2) - \mu_n^2} & \text{if } 0 \le n < P, \\ \gamma_n = -\sqrt{\mu_n^2 - (\omega^2/c^2)}, & \text{if } n \ge P, \end{array} \right. \qquad (28) $$
the modes are then given by the functions

$$v_n^{\pm}(x, z) = g_n(z)e^{\pm\gamma_n x}, \qquad (29)$$

where the ($g_n$) form a complete basis of $L^2(-H, 0)$. The modes $v_n^+$ are either propagating (for $n < P$) or evanescent (for $n \ge P$) in the right direction of $\Omega_0$, so that there are outgoing to the right. The modes $v_n^-$ are either propagating (for $n < P$) or evanescent (for $n \ge P$) in the left direction of $\Omega_0$, so that there are outgoing to the left.

The problem (23) in the unbounded domain $\Omega_0$ ($\theta = 0$) is equivalent to the following problem set in the bounded domain $\Omega_0^R$: find $u \in H^1(\Omega_0^R)$ such that

$$ \left\{ \begin{array}{l} \Delta u + (\omega^2/c^2)u = \begin{array}{c} 0 \\ u \end{array} \quad \text{in } \Omega_0^R, \\ \partial_z u = -i\omega\chi \quad \text{on } \Gamma_0^R, \\ \pm\partial_x u - T_{\pm}^a u = \begin{array}{c} 0 \\ 0 \end{array} \quad \text{on } \Sigma_{\pm R}, \end{array} \qquad (30) $$

where the Dirichlet-To-Neumann operators $T_{\pm}^a$ are defined as

$$ \left\{ \begin{array}{l} T_{\pm}^a : H^{1/2}(\Sigma_{\pm R}) \to \tilde{H}^{-1/2}(\Sigma_{\pm R}) \\ \varphi \mapsto \sum_{n \in \mathbb{N}} \gamma_n(\varphi, g_n) L^2(-H, 0)g_n, \end{array} \right. \qquad (31) $$

the $\gamma_n$ are given by (28) and the $g_n$ are defined by (27). Introducing a similar weak formulation as (13), we get

**Theorem 3.2.** The problem (30) has a unique solution.

Now let us consider the problem (4) when $g \to 0$. In this view, we introduce the corresponding fundamental solution in the waveguide $\Omega_0$, that is for some $N(x', z') \in \Omega_0$, the solution $G^a(\cdot, N) \in L_{\text{loc}}^2(\Omega_0)$ to the problem

$$ \left\{ \begin{array}{l} -\Delta G^a(\cdot, N) - (\omega^2/c^2)G^a(\cdot, N) = \delta_N \quad \text{in } \Omega_0, \\ G^a(\cdot, N) = 0 \quad \text{on } \Gamma_0, \\ \partial_z G^a(\cdot, N) = 0 \quad \text{on } \Gamma_{-H}, \\ G^a(\cdot, N) \text{ is outgoing.} \end{array} \right. $$

The expression of $G^a$ is given, for $M(x, z)$, by

$$ G^a(M, N) = - \sum_{n \in \mathbb{N}} \frac{1}{2\gamma_n} e^{\gamma_n |x-x'|} g_n(z) g_n(z'), \qquad (32) $$

where the $g_n$ and the $\gamma_n$ are given by (27) and (28), respectively. As in the case of the complete model, we obtain that the solution $u$ to the problem (4) in the particular case $g \to 0$ also satisfies in $\Omega_\theta^R$ the problem:

$$ \left\{ \begin{array}{l} \Delta u + (\omega^2/c^2)u = \begin{array}{c} 0 \\ u \end{array} \quad \text{in } \Omega_\theta^R, \\ \partial_\nu u = \begin{array}{c} 0 \\ 0 \end{array} \quad \text{on } \Gamma_0^R, \\ -\partial_x u = T_-^a u - 2\partial_x G^a(\cdot, N) \quad \text{on } \Sigma_{-R}, \\ \partial_x u = T_+^a u \quad \text{on } \Sigma_R. \end{array} \right. \qquad (33) $$
We have the following well-posedness result in the Fredholm sense.

**Theorem 3.3.** *Assume that the problem (33) has at most one solution. Then such problem is well-posed.*

**Remark 2.** Comparing the gravity modes (25) and the acoustic modes (28) allows us to anticipate the influence of the model on the resolution of the inverse problems, that is the tsunami identification and the bathymetry problem. Let us restrict to the right-going modes. In the gravity model we observe that the only one propagating mode $u_0^+$ is exponentially decaying to the bottom direction, while the evanescent modes $u_n^+$, $n \ge 1$, are all oscillating in the vertical direction. We hence expect that the information related to the bottom will be hard to identify from the knowledge of the free surface. Since $v_0$ is an increasing function with respect to $\omega$, we expect that the larger is $\omega$, the harder will be the identification. In the acoustic model, the modes $v_n^+$ are all oscillating in the vertical direction, notably the $P$ propagating ones. Having in mind that such number $P$ is increasing with respect to the frequency $\omega$, we expect that the inverse problems will be easier to solve in the acoustic case, in particular at high frequencies. Similarly, moving from the gravity case to the complete case, we observe in view of (9) that we increase the number of propagating modes, the first one being exponentially decreasing towards the bottom of the ocean, the other ones being oscillating in the depth direction. Of course, the resolution of the inverse problems will also benefit from this effect, in particular at high frequencies.

# 4 Setting of the inverse problems

## 4.1 The tsunami identification problem

Assume that for some unknown function $\chi \in L^2(\Gamma_{-H}^R)$, where $\chi$ is compactly supported in $(-R, R)$, the velocity potential $u$ is the solution to the problem (11), which is well-posed from Theorem 2.2. The inverse problem consists in retrieving the solicitation $\chi$ from the measurement of the induced free surface perturbation $\eta$. In this paper, we try to approach a concrete situation where the measured data correspond to a real and thus complicated model, while we solve the inverse problem from those data with the help of a simpler model. In this sense we avoid an inverse crime. More precisely, we will try to solve the inverse problem by using either of the three models: gravity, acoustic or complete. Let us first begin with the gravity model. The free surface perturbation is given by $\eta = (i\omega/g)u$ on $\Gamma_0^R$, such measurement being noisy in practice. It is readily seen that such linear inverse problem is ill-posed. Indeed, having in mind the Robin boundary condition $\partial_z u = (\omega^2/g)u$ on $\Gamma_0^R$, the inverse problem amounts to solve a Cauchy problem for the Laplace
equation in $\Omega_0^R$, the Cauchy data being $(u, \partial_z u)$ on $\Gamma_0^R$. Such problem is well-known to be severely ill-posed. In the case of the acoustic model, we also want to retrieve the solicitation $\chi$ from the measurement of the free surface perturbation $\eta$, but since $u$ satisfies a homogeneous Dirichlet condition on $\Gamma_0^R$, we now use the formula $\eta = -(1/i\omega)\partial_z u$ on $\Gamma_0^R$. Again we are brought back to a Cauchy problem for the Helmholtz equation, which is strongly ill-posed. If we use the complete model to solve the inverse problem, we use the relationship $\eta = (i\omega/g)u$ on $\Gamma_0^R$, like in the gravity model.

**Remark 3.** It is important to note that while the data $\eta$ of the inverse problem is proportional to the trace of $u$ on $\Gamma_0^R$ (Dirichlet data) in the complete/gravity models, it is proportional to the normal derivative $\partial_z u$ on $\Gamma_0^R$ (Neumann data) in the acoustic model. It will make a significant difference in the presence of noise as we will see hereafter.

Note that for the three models, there is at most one solution $\chi$ associated with the measurement $\eta$, this is a consequence of the uniqueness related to the Cauchy problem for the Laplace or Helmholtz equation. Here is a proof in the complete case. The other cases are similar.

**Theorem 4.1.** Let us consider two functions $\chi_1, \chi_2 \in L^2(\Gamma_0^R)$, and let us denote $u_j, j = 1, 2$, the corresponding solutions to the problem (11), as well as the corresponding free surface perturbations $\eta_j = (i\omega/g)u_j$ on $\Gamma_0^R$, $j = 1, 2$. If $\eta_1 = \eta_2$, then $\chi_1 = \chi_2$.

*Proof.* That $\eta_1 = \eta_2$ implies that $u_1|_{\Gamma_0^R} = u_2|_{\Gamma_0^R}$. Since for $j = 1, 2$, $u_j$ satisfies the boundary condition $\partial_z u_j - (\omega^2/g)u_j = 0$ on $\Gamma_0^R$ as well as the Helmholtz equation $\Delta u_j + (\omega^2/c^2)u_j = 0$ in $\Omega_0^R$, uniqueness for the Cauchy problem implies that $u_1 = u_2$ in $\Omega$. Since $\chi_j = (1/i\omega)\partial_z u_j$, we finally conclude that $\chi_1 = \chi_2$. $\square$

## 4.2 The bathymetry problem

The goal of bathymetry is to recover some unknown function $\theta \in C^0([-R, R])$ which is compactly supported in $(-R, R)$ from the free surface perturbation $\eta$ on $\Gamma_0^R$ induced by the point source $G^c(\cdot, N)$. Let us denote $u$ the solution to the problem (18) in the domain $\Omega_\theta^R$. Theorem 2.3 implies that such solution $u$ is well-defined in the absence of trapped solutions. If we try to solve the inverse problem by using the gravity or the complete model, we exploit the boundary condition $\eta = (i\omega/g)u$ on $\Gamma_0^R$. If we use the acoustic model, the goal is the same but we use the formula $\eta = -(1/i\omega)\partial_z u$ on $\Gamma_0^R$. As a geometric inverse problem, bathymetry is a non-linear inverse problem and is of course ill-posed. Uniqueness of the bathymetry problem follows the same lines as the classical Schiffer proof (see for example [15]). For the sake of self-containment, we offer a proof in the gravity case in the case of a smooth function $\theta$. The other cases are similar.
**Theorem 4.2.** Let us consider two $C^3$ functions $\theta_1, \theta_2$ which are compactly supported in $(-R, R)$, and let us denote $u_j, j = 1, 2$, the corresponding solutions to the problem (18) with $c \to +\infty$, as well as the corresponding free surface perturbations $\eta_j = (i\omega/g)u_j$ on $\Gamma_0^R, j = 1, 2$. We assume that $\eta_1$ and $\eta_2$ are not constant functions. If $\eta_1 = \eta_2$, then $\theta_1 = \theta_2$.

*Proof.* Since $u_1$ and $u_2$ have the same trace on $\Gamma_0^R$ and satisfy the same Robin boundary condition on $\Gamma_0^R$, we have $(u_1|_{\Gamma_0^R}, \partial_z u_1|_{\Gamma_0^R}) = (u_2|_{\Gamma_0^R}, \partial_z u_2|_{\Gamma_0^R})$. From uniqueness for the Cauchy problem for the Laplace equation, we get that $u_1 = u_2$ in $I_{12}^R := \Omega_{\theta_1}^R \cap \Omega_{\theta_2}^R$. Assume that $\Omega_{\theta_2} \not\subset \Omega_{\theta_1}$ and let us consider $\mathcal{R} := U_{12}^R \setminus \Omega_{\theta_1}^R$, where $U_{12}^R := \Omega_{\theta_1}^R \cup \Omega_{\theta_2}^R$, which is not empty. We have $\partial \mathcal{R} \subset (\partial \Omega_{\theta_1}^R \cap \partial I_{12}^R) \cup \partial \Omega_{\theta_2}^R$. By denoting $\nu$ the outward unit normal to $\mathcal{R}$, on $\partial \Omega_{\theta_1}^R \cap \partial I_{12}^R$, we have $\partial_\nu u_1 = \partial_\nu u_2 = 0$, while on $\partial \Omega_{\theta_2}^R$, we have $\partial_\nu u_2 = 0$. As a result, the function $u_2$ satisfies $\Delta u_2 = 0$ in $\mathcal{R}$ with the boundary condition $\partial_\nu u_2 = 0$ on $\partial \mathcal{R} \cap \Gamma_{-H+\theta_1}^R$ and $\partial_\nu u_2 = 0$ on $\partial \mathcal{R} \cap \Gamma_{-H+\theta_2}^R$. We observe here that the domain $\mathcal{R}$ is not necessarily Lipschitz, which means that the standard Green formula in $\mathcal{R}$ does not hold. To cope with this problem we follow the argument of [2]. Since for $j = 1, 2$ the function $\theta_j$ is in $C^3([-R, R])$, the corresponding solution $u_j$ is $C^2$ in a volumic neighborhood of a subset of the bottom $\Gamma_{-H+\theta_j}^R$ which does not meet the sections $\Sigma_{\pm R}$ (see [23]). Hence $u_2 \in C^2(\bar{\mathcal{R}})$. The regularity of $u_2$ and that the perimeter of $\mathcal{R}$ is finite imply the generalized Green formula (see [21]):

$$ \int_{\mathcal{R}} \Delta u_2 \overline{u_2} dx = - \int_{\mathcal{R}} |\nabla u_2|^2 dx + \int_{\partial \mathcal{R}} \partial_\nu u_2 \overline{u_2} ds, $$

where $\partial^* \mathcal{R}$ is the reduced boundary in the sense of De Georgi. We get that $\nabla u_2 = 0$ in $\mathcal{R}$, that is $u_2 = C$ in $\mathcal{R}$, where $C$ is a constant. From unique continuation applied to the function $u_2 - C$, we get that $u_2 = C$ in $\Omega_{\theta_2}^R$, which implies that $\eta_2$ is constant on $\Gamma_0^R$. This is a contradiction. $\square$

# 5 An abstract framework to solve ill-posed Cauchy problems

In this section we introduce an abstract framework which will be well-adapted to the regularization of the two inverse problems presented above (tsunami identification and bathymetry) and for the three different models (gravity, acoustic and complete models). This abstract framework is already presented in [11] in the case of real Hilbert spaces, and we have to adapt it to the case of complex Hilbert spaces. We then consider three complex Hilbert spaces $V, M$ and $O$, equipped with the scalar products $(\cdot, \cdot)_V, (\cdot, \cdot)_M$ and $(\cdot, \cdot)_O$, respectively. We also consider two continuous operators $B: V \to M$ and $C: V \to O$, as well as $A: V \to M \times O$ defined with the help of $B$ and
C by $Av = (Bv, Cv)$, for all $v \in V$. We assume that A is injective and has a dense range. However, it is not supposed to be onto, which implies that the problem: for $F = (F_1, F_2) \in M \times O$, find $u \in V$ such that $Au = F$ may have no solution. For $\varepsilon > 0$, the standard Tikhonov regularization consists in the well-posed weak formulation: find $u_\varepsilon \in V$ such that

$$ (Au_\varepsilon, Av)_{M \times O} + \varepsilon(u_\varepsilon, v)_V = (F, Av)_{M \times O}, \quad \forall v \in V. \quad (34) $$

The idea of the mixed formulation of the Tikhonov regularization consists in introducing an auxiliary unknown $\lambda_\varepsilon \in M$.

**Proposition 1.** For $F = (F_1, F_2) \in M \times O$, $u_\varepsilon \in V$ is the solution to the problem (34) if and only if $(u_\varepsilon, \lambda_\varepsilon := Bu_\varepsilon - F_1) \in V \times M$ is the solution to the following weak formulation:

$$ \left\{ \begin{array}{l} \varepsilon(u_\varepsilon, v)_V + \overline{b(v, \lambda_\varepsilon)} + (Cu_\varepsilon, Cv)_O = (F_2, Cv)_O, \quad \forall v \in V, \\ b(u_\varepsilon, \mu) - (\lambda_\varepsilon, \mu)_M = (F_1, \mu)_M, \quad \forall \mu \in M, \end{array} \right. \quad (35) $$

where the sesquilinear form $b$ on $V \times M$ is defined by

$$ b(v, \mu) = (Bv, \mu)_M, \quad \forall (v, \mu) \in V \times M. \quad (36) $$

Such weak formulation is well-posed for all $\varepsilon > 0$ and if there exists $u \in V$ such that $Au = F$, then $u_\varepsilon \to u$ in $V$ and $\lambda_\varepsilon \to 0$ in $M$ when $\varepsilon \to 0$.

*Proof.* That $u_\varepsilon$ satisfies (34) is equivalent to the fact that

$$ (Bu_\varepsilon, Bv)_M + (Cu_\varepsilon, Cv)_O + \varepsilon(u_\varepsilon, v)_V = (F_1, Bv)_M + (F_2, Cv), \quad \forall v \in V. $$

Introducing $\lambda_\varepsilon := Bu_\varepsilon - F_1$, we get that the above weak formulation is equivalent to find $(u_\varepsilon, \lambda_\varepsilon) \in V \times M$ such that

$$ \left\{ \begin{array}{l} \varepsilon(u_\varepsilon, v)_V + \overline{(Bv, \lambda_\varepsilon)_M} + (Cu_\varepsilon, Cv)_O = (F_2, Cv)_O, \quad \forall v \in V, \\ (Bu_\varepsilon, \mu)_M - (\lambda_\varepsilon, \mu)_M = (F_1, \mu)_M, \quad \forall \mu \in M, \end{array} \right. $$

that is (35). Since the problem (34) is well-posed, this is also true for the problem (35). Besides, if there exists $u \in V$ such that $Au = F$, a well-known property of the Tikhonov regularization is that $u_\varepsilon \to u$ in $V$ when $\varepsilon \to 0$. Then $\lambda_\varepsilon = Bu_\varepsilon - F_1 \to Bu - F_1 = 0$ when $\varepsilon \to 0$. $\square$

In the inverse problems described in the previous section, the data $F = (F_1, F_2)$ is noisy in practice, so that we measure some data $F^\delta = (F_1^\delta, F_2^\delta) \in M \times O$, which satisfies the classical assumption

$$ \|F^\delta - F\|_{M \times O} \le \delta < \|F^\delta\|_{M \times O}. \quad (37) $$

In the applications hereafter, we will see that only one component among $F_1$ or $F_2$ is noisy, the other component being uncontaminated. In the first case

$$ \|F^\delta - F\|_{M \times O} = \|F_1^\delta - F_1\|_M $$
while in the second case

$$ \|F^{\delta} - F\|_{M \times O} = \|F_2^{\delta} - F_2\|_{O}. $$

The Morozov principle offers a rigorously defined way of choosing the regularization parameter $\varepsilon$ in the Tikhonov formulation (34), or equivalently in the mixed Tikhonov regularization (35), as a function of the amplitude of noise $\delta$. We have the following result.

**Theorem 5.1.** Let us consider some data $F^{\delta} \in M \times O$ satisfying (37), and $u_{\varepsilon}^{\delta} \in V$ the solution to the problem (34) with data $F$ replaced by $F^{\delta}$. There exists a unique $\varepsilon > 0$ such that

$$ \|Au_{\varepsilon}^{\delta} - F^{\delta}\|_{M \times O} = \delta. \qquad (38) $$

The proof of Theorem 5.1 is omitted since it is a slight adaptation to complex Hilbert spaces of Theorem 2.1 in [11]. We note that the equality (38) is equivalent to

$$ \|\lambda_{\varepsilon}^{\delta}\|_{M}^{2} + \|Cu_{\varepsilon}^{\delta} - F_{2}^{\delta}\|_{O}^{2} = \delta^{2}, $$

where $(u_{\varepsilon}^{\delta}, \lambda_{\varepsilon}^{\delta}) \in V \times M$ is the solution to the problem (35) for data $F^{\delta} = (F_{1}^{\delta}, F_{2}^{\delta}) \in M \times O$ instead of $F = (F_{1}, F_{2})$.

The solution $u_{\varepsilon}^{\delta} \in V$ to the problem (34) with data $F^{\delta}$ satisfies

$$ A^{*}Au_{\varepsilon}^{\delta} + \varepsilon u_{\varepsilon}^{\delta} = A^{*}F^{\delta}, \qquad (39) $$

where $A^{*}: M \times O \to V$ is the adjoint operator of $A$. It remains to find a strategy to compute the value $\varepsilon(\delta)$ defined by Theorem 5.1, as well as the corresponding solution $u^{\delta} := u_{\varepsilon(\delta)}^{\delta}$. One way to do that is to solve the following minimization problem in $H := M \times O$:

$$ \inf_{q \in H} J^{\delta}(q), \quad J^{\delta}(q) := \frac{1}{2} \|A^{*}q\|_{V}^{2} + \delta \|q\|_{H} - \mathcal{R}e \{(F^{\delta}, q)_{H}\}. \qquad (40) $$

We have the following proposition.

**Proposition 2.** The problem (40) has a unique solution $p^{\delta} \in H$. In addition, we have

$$ \varepsilon(\delta) = \frac{\delta}{\|p^{\delta}\|_{H}}, \quad u^{\delta} = A^{*}p^{\delta}. $$

*Proof.* Let us prove that the functional $J^{\delta}$ is coercive, that is $J^{\delta}(q) \to +\infty$ when $\|q\|_{H} \to +\infty$. Assume on the contrary that there exists some constant $C > 0$ and a sequence $(q_n)_{n \in \mathbb{N}}$ of elements in $H$ such that $\alpha_n = \|q_n\|_{H} \to +\infty$ while $J^{\delta}(q_n) \le C$. Let us define $z_n = q_n/\|q_n\|_{H}$. Since $\|z_n\|_{H} = 1$, there exists a subsequence of $(z_n)_n$, still denoted $(z_n)_n$, such that $z_n \to z$ in $H$. We have

$$ \frac{\alpha_n^2}{2} \|A^*z_n\|_V^2 + \alpha_n \delta - \alpha_n \mathcal{R}e \{(F^\delta, z_n)_H\} \le C, \qquad (41) $$
hence

$$ \frac{1}{2} \|A^*z_n\|_V^2 \le \frac{\|F^\delta\|_H}{\alpha_n} + \frac{C}{\alpha_n^2}, $$

which implies that $A^*z_n \to 0$ in $V$ and since $A^*z_n \to A^*z$ in $V$, we get that $A^*z = 0$. That $A$ has a dense range implies that $A^*$ is injective, hence $z = 0$, that is $z_n \to 0$ in $H$. Another consequence of (41) is that

$$ C \ge \alpha_n(\delta - \mathcal{R}e\{(F^\delta, z_n)_H\}) \to +\infty, $$

which is a contradiction. The functional $J^\delta$ is coercive, continuous and strictly convex, which implies that the problem (40) has a unique solution $p^\delta$. Let us check that $p^\delta \ne 0$. Actually, let us take $z^\delta = F^\delta/\|F^\delta\|_H$. For $\varepsilon > 0$, we have

$$ J^\delta(\varepsilon z^\delta) = \frac{\varepsilon^2}{2} \|A^* z^\delta\|_V^2 + \varepsilon (\delta - \|F^\delta\|_H). $$

For small $\varepsilon$, the real number $J^\delta(\varepsilon z^\delta)$ has the sign of $\delta - \|F^\delta\|_H < 0$, hence there exists $q \in H$ such that $J^\delta(q) < 0 = J^\delta(0)$, which implies that $p^\delta$ is not zero. The functional $J^\delta$ is differentiable in the sense of Fréchet at $p \ne 0$ and its derivative is given by

$$ \langle DJ^\delta(p), q \rangle_{H', H} = \mathcal{R}e \left\{ (AA^*p, q)_H + \frac{\delta}{\|p\|_H} (p, q)_H - (F^\delta, q)_H \right\}, \quad \forall q \in H. $$

Writing that $\langle DJ^\delta(p^\delta), q \rangle_{H', H} = 0$ for any $q$ and $iq$ in $H$, we finally get that

$$ AA^*p^\delta + \frac{\delta}{\|p^\delta\|_H} p^\delta = F^\delta, $$

Denoting $w^\delta := A^*p^\delta \in V$, we obtain that

$$ A^*Aw^\delta + \frac{\delta}{\|p^\delta\|_H} w^\delta = A^*F^\delta, \quad \|Aw^\delta - F^\delta\|_H = \delta. $$

This means that for $\varepsilon = \delta/\|p^\delta\|_H$, $w^\delta$ satisfies both (39) and (38), that is $w^\delta = u^\delta$. In other words, such value of $\varepsilon$ corresponds to the Morozov choice in the Tikhonov regularization. $\square$

We observe that the functional $J^\delta$ given by (40) is convex but is not quadratic because of the second term $\delta\|q\|_H$. In practice, the minimization of $J^\delta$ by a gradient method is likely to be more difficult than if it were quadratic. In [16], the author proposed a strategy to circumvent this issue. It consists in the following iterative algorithm, where the solution of the initial non-quadratic minimization problem (40) is approached by a sequence formed by the minimizers of well-defined quadratic functionals.
**Algorithm 5.2.** 1. Let $\varepsilon_n > 0$ be given

2. Compute the unique solution $p_n$ of the minimization problem

$$ \inf_{q \in H} L_n(q), \quad L_n(q) := \frac{1}{2} \|A^*q\|_V^2 + \frac{1}{2} \varepsilon_n \|q\|_H^2 - \mathcal{R}e \left\{ (F^\delta, q)_H \right\} \quad (42) $$

3. Compute $\varepsilon_{n+1} = \delta/\|p_n\|_H$ and go to step 1 for $n \to n+1$

Comparing (40) and (42), the second non-quadratic term $\delta\|q\|_H$ was replaced by the quadratic term $\varepsilon_n\|q\|_H^2/2$. In [16], the following proposition in proved.

**Proposition 3.** Let $(p_n)_{n \in \mathbb{N}}$ be the sequence of elements of $H$ defined by the algorithm 5.2. Then $p_n \to p^\delta$ in $H$ when $n \to +\infty$ and $u_n := A^*p_n \to u^\delta$ in $V$ when $n \to +\infty$.

We now establish that the minimization problem (42) is equivalent to a weak formulation satisfied by $u_n \in V$ and $p_n = (\lambda_n, h_n) \in M \times O$ for all $n \in \mathbb{N}$.

**Proposition 4.** The triple $(u_n, \lambda_n, h_n) \in V \times M \times O$ is given by the well-posed weak formulation: for all $(v, \mu, h) \in V \times M \times O$,

$$ \left\{ \begin{array}{l} (u_n, v)_V - \overline{b(v, \lambda_n)} - (h_n, Cv)_O = 0, \quad \forall v \in V, \\ b(u_n, \mu) + \varepsilon_n(\lambda_n, \mu)_M = (F_1^\delta, \mu)_M, \quad \forall \mu \in M, \\ (Cu_n, h)_O + \varepsilon_n(h_n, h)_O = (F_2^\delta, h)_O, \quad \forall h \in O. \end{array} \right. \quad (43) $$

*Proof.* Let us first notice that $A^*(\lambda, h) = A^*p = B^*\lambda + C^*h$ for all $p = (\lambda, h) \in M \times O$, hence,

$$ L_n(\lambda, h) = \frac{1}{2} \|B^*\lambda + C^*h\|_V^2 + \frac{\varepsilon_n}{2} (\|\lambda\|_M^2 + \|h\|_O^2) - \mathcal{R}e \left\{ (F_1^\delta, \lambda)_M + (F_2^\delta, h)_O \right\}. $$

The minimizer $(\lambda_n, h_n)$ of $L_n$ satisfies

$$ \left\{ \begin{array}{l} (B^*\lambda_n + C^*h_n, B^*\mu)_V + \varepsilon_n(\lambda_n, \mu)_M - (F_1^\delta, \mu)_M = 0, \quad \forall \mu \in M, \\ (B^*\lambda_n + C^*h_n, C^*h)_V + \varepsilon_n(h_n, h)_O - (F_2^\delta, h)_O = 0, \quad \forall h \in O. \end{array} \right. $$

Since $u_n = B^*\lambda_n + C^*h_n$, we get that

$$ \left\{ \begin{array}{l} (u_n, v)_V - (B^*\lambda_n, v)_V - (C^*h_n, v)_V = 0, \quad \forall v \in V, \\ (u_n, B^*\mu)_V + \varepsilon_n(\lambda_n, \mu)_M - (F_1^\delta, \mu)_M = 0, \quad \forall \mu \in M, \\ (u_n, C^*h)_V + \varepsilon_n(h_n, h)_O - (F_2^\delta, h)_O = 0, \quad \forall h \in O. \end{array} \right. $$

By using the definition of the operators $B^*$, $C^*$ and the sesquilinear $b$, we obtain the system (43), which completes the proof of the proposition. $\square$
# 6 Application to the tsunami identification

The goal of the section is to apply the abstract setting of the previous section to the tsunami identification problem, both for the gravity and the acoustic models. The case of the complete model is omitted, since it is very close to the gravity one.

## 6.1 The gravity model

Let us start with the gravity model. The potential $u$ satisfies the problem

$$ egin{cases} 
\Delta u = 0 & \text{in } \Omega_0^R, \\
\partial_z u - (\omega^2/g)u = 0 & \text{on } \Gamma_0^R, \\
\pm \partial_x u - T_{\pm}^g u = 0 & \text{on } \Sigma_{\pm R},
\end{cases} 
\quad (44) $$

with additionally

$$ u = F_2 := (g/i\omega)\eta \quad \text{on } \Gamma_0^R, \quad (45) $$

where $\eta$ is the measured free surface perturbation. In (44), we recall that the operators $T_{\pm}^g$ coincide with the operators $T_{\pm}^c$ defined by (12) when $c \to +\infty$, that is the functions $f_n$ and the numbers $\beta_n$ are defined by (24) and (25), respectively. The goal is to find an approximation of $u$ in $\Omega_0^R$ from the data $F_2$, which is noisy in general, using the strategy described in the previous section. We specify the Hilbert spaces $V$, $M$ and $O$ as

$$ V = H^1(\Omega_0^R), \quad M = \{\mu \in H^1(\Omega_0^R), \mu|_{\Gamma_{-H}^R} = 0\}, \quad O = L^2(\Gamma_0^R). \quad (46) $$

The norm of space $M$ is chosen as the $H^1(\Omega_0^R)$ semi-norm, thanks to the Poincaré inequality. We also introduce the continuous sesquilinear form on $V \times M$:

$$ b(v, \mu) = \int_{\Omega_0^R} \nabla v \cdot \nabla \bar{\mu} \,dx dz - (\omega^2/g) \int_{\Gamma_0^R} v \bar{\mu} \,dx \\ - \langle T_+^g v, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} - \langle T_-^g v, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \quad (47) $$

and the corresponding operator $B: V \to M$ defined by (36) from such sesquilinear form. The operator $C: V \to O$ is chosen as the trace operator. This enables us to define the operator $A: V \to M \times O$ as $A = (B, C)$. By construction $F_1 = 0$, so that $F = (0, F_2) \in M \times O$. The following lemma establishes a connection between the tsunami identification and the abstract framework of the previous section.

**Lemma 6.1.** The potential $u \in H^1(\Omega_0^R)$ satisfies the problem (44) (45) if and only if $u$ satisfies $Au = F$.
Proof. Assume that $u \in H^1(\Omega_0^R)$ satisfies $Au = F$. We hence have $Bu = 0$ and $Cu = F_2$. That $Cu = F_2$ is equivalent to $u|_{\Gamma_0^R} = F_2$. That $Bu = 0$ is equivalent to $b(u, \mu) = 0$ for all $\mu \in M$, that is

$$ \int_{\Omega_0^R} \nabla u \cdot \nabla \bar{\mu} \,dx\,dz - (\omega^2/g) \int_{\Gamma_0^R} u \bar{\mu} \,dx - \langle T_+^g u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\ - \langle T_-^g u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} = 0, \quad \forall \mu \in M. $$

Let us first choose $\mu = \varphi \in C_0^\infty(\Omega_0^R)$. We obtain that $\Delta u = 0$ in $\Omega_0^R$. By using the Green formula, we get

$$ \int_{\Omega_0^R} \nabla u \cdot \nabla \bar{\mu} \,dx\,dz = - \int_{\Omega_0^R} \Delta u \bar{\mu} \,dx\,dz + \langle \partial_\nu u, \mu \rangle_{H^{-1/2}(\partial \Omega_0^R), H^{1/2}(\partial \Omega_0^R)}, \quad \forall \mu \in M, $$
whence

$$ \langle \partial_{\nu} u, \bar{\mu} \rangle_{H^{-1/2}(\partial \Omega_0^R), H^{1/2}(\partial \Omega_0^R)} - (\omega^2/g) \int_{\Gamma_0^R} u \bar{\mu} \,dx - \langle T_+^g u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\ - \langle T_-^g u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} = 0, \quad \forall \mu \in M. $$

This yields $\partial_z u = (\omega^2/g)u$ on $\Gamma_0^R$ and that $\pm \partial_x u = T_\pm^g u$ on $\Sigma_{\pm R}$. We eventually establish that $u$ satisfies (44) and (45). The converse statement is obtained analogously. $\square$

To apply the results of the previous section, we have to check that the operator $A$ is injective and has a dense range.

**Lemma 6.2.** *The operator A is injective and has a dense range.*

Proof. Let us consider $u \in H^1(\Omega_0^R)$ such that $Au = 0$. From the proof of Lemma 6.1, we have that $\Delta u = 0$ in $\Omega_0^R$, $u = 0$ on $\Gamma_0^R$ and $\partial_z u = (\omega^2/g)u$ on $\Gamma_0^R$, that is $(u, \partial_z u) = (0, 0)$ on $\Gamma_0^R$. From uniqueness of the Cauchy problem for the Laplace equation, we have $u = 0$ in $\Omega_0^R$, which shows the injectivity of the operator $A$.

Now let us assume that $(\mu, h) \in M \times L^2(\Gamma_0^R)$ satisfies

$$ (Av, (\mu, h))_{M \times L^2(\Gamma_0^R)} = 0, \quad \forall v \in H^1(\Omega_0^R). $$

This implies that

$$ \int_{\Omega_0^R} \nabla v \cdot \nabla \bar{\mu} \,dx\,dz - (\omega^2/g) \int_{\Gamma_0^R} v \bar{\mu} \,dx - \langle T_+^g v, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\ - \langle T_-^g v, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} + \int_{\Gamma_0^R} v \bar{h} \,dx = 0, \quad \forall v \in H^1(\Omega_0^R). $$

Choosing first $v = \varphi \in C_0^\infty(\Omega_0^R)$, we get $\Delta\mu = 0$ in $\Omega_0^R$. Using the Green formula, we obtain that $\partial_z\mu = 0$ on $\Gamma_{-H}^R$ and $\partial_z\mu - (\omega^2/g)\mu + h = 0$ on $\Gamma_0^R$. Since we also have $\mu = 0$ on $\Gamma_{-H}^R$, uniqueness of the Cauchy problem for the Laplace equation implies that $\mu = 0$ in $\Omega_0^R$, which in turn implies that $h = 0$ in $L^2(\Gamma_0^R)$. We conclude that the operator $A$ has a dense range. $\square$
We can then apply all the results of the abstract framework described in the previous section to that specific example. In particular, the weak formulation (35) becomes (note that $F_1 = 0$): find $(u_\varepsilon, \lambda_\varepsilon) \in H^1(\Omega_0^R) \times M$, with $M = \{\mu \in H^1(\Omega_0^R), \mu|_{\Gamma_{-H}^R} = 0\}$, such that

$$ egin{cases} 
\varepsilon(u_\varepsilon, v)_{H^1(\Omega_0^R)} + \int_{\Omega_0^R} \nabla\lambda_\varepsilon \cdot \nabla\bar{v} \,dx\,dz - (\omega^2/g) \int_{\Gamma_0^R} \lambda_\varepsilon \bar{v} \,dx \\ 
-\langle \overline{T_+^g} \lambda_\varepsilon, \bar{v} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} - \langle \overline{T_-^g} \lambda_\varepsilon, \bar{v} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \\ 
+ \int_{\Gamma_0^R} u_\varepsilon \bar{v} \,dx = \int_{\Gamma_0^R} F_2 \bar{v} \,dx, \quad \forall v \in H^1(\Omega_0^R), \\ 
\int_{\Omega_0^R} \nabla u_\varepsilon \cdot \nabla\bar{\mu} \,dx\,dz - (\omega^2/g) \int_{\Gamma_0^R} u_\varepsilon \bar{\mu} \,dx \\ 
-\langle \overline{T_+^g} u_\varepsilon, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} - \langle \overline{T_-^g} u_\varepsilon, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \\ 
- \int_{\Omega_0^R} \nabla\lambda_\varepsilon \cdot \nabla\bar{\mu} \,dx\,dz = 0, \quad \forall \mu \in M. 
\end{cases} 
\quad (48) $$

In the presence of a noisy data $F_2^\delta$ such that

$$ \|F_2^\delta - F_2\|_{L^2(\Gamma_0^R)} \le \delta < \|F_2^\delta\|_{L^2(\Gamma_0^R)}, \quad (49) $$

which, since $F_1 = 0$, amounts to (37), we apply the algorithm 5.2 to obtain the Morozov solution $u^\delta \in H^1(\Omega_0^R)$. This requires to solve for any $n \in \mathbb{N}$ the problem (43), which becomes here: find $(u_n, \lambda_n, h_n) \in H^1(\Omega_0^R) \times M \times L^2(\Gamma_0^R)$ such that

$$ egin{cases} 
(u_n, v)_{H^1(\Omega_0^R)} - \int_{\Omega_0^R} \nabla \lambda_n \cdot \nabla \bar{v} \,dx\,dz + (\omega^2/g) \int_{\Gamma_0^R} \lambda_n \bar{v} \,dx \\ 
+ \langle \overline{T_+^g} \lambda_n, \bar{v} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} + \langle \overline{T_-^g} \lambda_n, \bar{v} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \\ 
- \int_{\Gamma_0^R} h_n \bar{v} \,dx = 0, \quad \forall v \in H^1(\Omega_0^R), \\ 
\int_{\Omega_0^R} \nabla u_n \cdot \nabla \bar{\mu} \,dx\,dz - (\omega^2/g) \int_{\Gamma_0^R} u_n \bar{\mu} \,dx \\ 
- \langle \overline{T_+^g} u_n, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} - \langle \overline{T_-^g} u_n, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \\ 
+ \varepsilon_n \int_{\Omega_0^R} \nabla \lambda_n \cdot \nabla \bar{\mu} \,dx\,dz = 0, \quad \forall \mu \in M, \\ 
\int_{\Gamma_0^R} u_n \bar{h} \,dx + \varepsilon_n \int_{\Gamma_0^R} h_n \bar{h} \,dx = \int_{\Gamma_0^R} F_2^\delta \bar{h} \,dx, \quad \forall h \in L^2(\Gamma_0^R). 
\end{cases} 
\quad (50) $$

The method that we propose to approximate the tsunami shape $\chi$ on $\Gamma_{-H}^R$ is the following: starting with an initial value $\varepsilon_0 > 0$, we successively solve
(50) and compute for $n \in \mathbb{N}$:

$$ \varepsilon_{n+1} = \frac{\delta}{\|p_n\|_H}, \quad \text{with} \quad \|p_n\|_H = \sqrt{\int_{\Omega_0^R} |\nabla \lambda_n|^2 dx dz + \int_{\Gamma_0^R} |h_n|^2 dx}. $$

For sufficiently large $n$, $u^\delta$ is approximated by $u_n$, which we still denote $u^\delta$ for simplicity. Lastly, in view of (11), the function $\chi$ is approximated by

$$ \chi^\delta = \frac{i}{\omega} \partial_z u^\delta. $$

## 6.2 The acoustic model

In the case of the acoustic model, the potential $u$ satisfies the problem

$$ \left\{ \begin{array}{lll} \Delta u + (\omega^2/c^2)u = 0 & \text{in } \Omega_0^R, \\ u = 0 & \text{on } \Gamma_0^R, \\ \pm \partial_x u - T_{\pm}^a u = 0 & \text{on } \Sigma_{\pm R}, \end{array} \right. \quad (51) $$

with additionally

$$ \partial_z u = f_1 := -i\omega\eta \quad \text{on } \Gamma_0^R, \quad (52) $$

where $\eta$ is the measured free surface perturbation. In (51), the operators $T_{\pm}^a$ are defined by (31). Again, the goal is to find an approximation of $u$ in $\Omega_0^R$ from the data $f_1$, which is noisy, by using the strategy described in the previous section. We specify the Hilbert spaces $V$, $M$ and $O$ as in (46). While the operator $C: H^1(\Omega_0^R) \to L^2(\Gamma_0^R)$ is the trace operator again, we introduce a new continuous sesquilinear form on $V \times M$:

$$ \begin{aligned} b(v, \mu) = & \int_{\Omega_0^R} \nabla v \cdot \nabla \bar{\mu} \,dx dz - (\omega^2/c^2) \int_{\Omega_0^R} v \bar{\mu} \,dx dz \\ & - \langle T_+^a u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} - \langle T_-^a u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} \end{aligned} $$

and the corresponding operator $B: V \to M$ defined by (36). We again form the operator $A: V \to M \times O$ as $A = (B, C)$. Let us now denote $L_1 \in M$ the unique element of $M$ which satisfies

$$ (L_1, \mu)_M = \int_{\Gamma_0^R} f_1 \bar{\mu} \,dx, \quad \forall \mu \in M. \quad (53) $$

We choose $F_1 = L_1$ and $F_2 = 0$, so that $F = (L_1, 0) \in M \times O$. The following lemma establishes a connection between the tsunami identification and the abstract framework of the previous section.

**Lemma 6.3.** The potential $u \in H^1(\Omega_0^R)$ satisfies the problem (51) (52) if and only if $u$ satisfies $Au = F$.
*Proof.* Let us consider $u \in H^1(\Omega_0^R)$ such that $Au = F$, that is $Bu = F_1$ and $Cu = 0$. That $Cu = 0$ is equivalent to $u = 0$ on $\Gamma_0^R$. That $Bu = F_1$ is equivalent to

$$ \int_{\Omega_0^R} \nabla u \cdot \nabla \bar{\mu} \, dx dz - (\omega^2/c^2) \int_{\Omega_0^R} u \bar{\mu} \, dx dz - \langle T_+^a u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_R), \tilde{H}^{1/2}(\Sigma_R)} \\ - \langle T_-^a u, \bar{\mu} \rangle_{H^{-1/2}(\Sigma_{-R}), \tilde{H}^{1/2}(\Sigma_{-R})} = \int_{\Gamma_0^R} f_1 \, \mu \, dx, \quad \forall \mu \in M. $$

Choosing first $\mu = \varphi \in C_0^\infty(\Omega_0^R)$ yields $\Delta u + (\omega^2/c^2)u = 0$ in $\Omega_0^R$. Then by using the Green formula, we obtain $\partial_z u = f_1$ on $\Gamma_0^R$ and $\pm \partial_x u = T_\pm^a u$ on $\Sigma_{\pm R}$, that is $u$ solves the problem (51) (52). The converse statement is obtained the same way. $\square$

Lemma 6.2 is also satisfied for our new specific operator $A$, so that we can again apply all the results of the abstract framework described in the previous section. The explicit forms of the systems (35) and (43) in the case of acoustics, which are not given here for brevity, are quite close to (48) and (50), respectively.

# 7 Application to bathymetry

We now apply the abstract setting of section 5 to the bathymetry problem, both for the gravity and the acoustic models. The case of the complete model is omitted, since it is very close to the gravity one. Contrary to the tsunami identification problem, we have to tackle a (non linear) geometric inverse problem. If we assume that the bottom perturbation $\theta$ is small, we propose to solve the bathymetry problem by approximating the boundary condition $\partial_\nu u = 0$ on the perturbed boundary $\Gamma_{-H+\theta}^R$ in (18) and in (33) by a generalized impedance boundary condition on the flat boundary $\Gamma_{-H}^R$, as proposed for example in [3] and [7] (see also the introduction of [9]).

## 7.1 The gravity model

Let us start with the gravity model. By reproducing the computations of [3] or [9], we can show that the first order approximation on $\Gamma_{-H}^R$ is $\partial_z u = 0$ while the second order approximation is

$$ \partial_z u = \partial_x (\theta \partial_x u) \quad \text{on} \quad \Gamma_{-H}^R, \qquad (54) $$

where $\|\theta\|_{L^\infty(-R,R)}$ is supposed to be small compared to the wavelength of the wave phenomenon at stake, that is $\lambda = 2\pi\sqrt{g\tilde{H}/\omega}$ in the gravity model.
In view of (18) with $c \to +\infty$, the velocity potential satisfies

$$ \left\{ \begin{array}{rcll} \Delta u & = & 0 & \text{in } \Omega_0^R, \\ \partial_z u - (\omega^2/g)u & = & 0 & \text{on } \Gamma_0^R, \\ -\partial_x u & = & T_-^g u - 2\partial_x \mathcal{G}^g(\cdot, N) & \text{on } \Sigma_{-R}, \\ \partial_x u & = & T_+^g u & \text{on } \Sigma_R \end{array} \right. \qquad (55) $$

and (45), where the operators $T_{\pm}^g$ are defined by (12) for $c \to +\infty$ and the point source $\mathcal{G}^g(\cdot, N)$ is given by (16) for $c \to +\infty$. We reuse the Hilbert spaces $V$, $M$ and $O$ as well as the operators $A$, $B$ and $C$ specified in section 6.1. We also define $G^g \in M$ such that

$$ (G^g, \mu)_M = -2 \int_{\Sigma_{-R}} \partial_x \mathcal{G}^g(\cdot, N) \bar{\mu} \, dy, \quad \forall \mu \in M. \qquad (56) $$

We choose $F_1 = G^g$ and $F = (G^g, F_2) \in M \times O$, with $G^g$ and $F_2$ given by (56) and (45), respectively. We easily check the following result.

**Lemma 7.1.** The function $u \in H^1(\Omega_0^R)$ satisfies the problem (55) (45) if and only if $u$ satisfies $Au = F$.

Such Lemma enables us to directly apply all the results of section 5. In particular, when the noisy data $F_2^\delta \in O$ is such that the $F^\delta = (G^g, F_2^\delta)$ satisfies (37), we again apply the algorithm 5.2 to obtain the Morozov solution $u^\delta \in H^1(\Omega_0^R)$.

In order to approximate $\theta$ by some function $\theta^\delta$ in the presence of noisy data, let us set

$$ \alpha := u^\delta|_{\Gamma_{-H}^R}, \quad \alpha_x := d\alpha/dx, \quad \alpha_z := \partial_z u^\delta|_{\Gamma_{-H}^R}, \qquad (57) $$

where the function $u^\delta$ has been obtained in the first step, and let us compute $\theta^\delta$ from the equation (54). From the strong equation satisfied by $u^\delta \in H^1(\Omega_0^R)$, we have $\Delta u^\delta \in L^2(\Omega_0^R)$, so that $\alpha \in H^{1/2}(\Gamma_{-H}^R)$, $\alpha_x \in H^{-1/2}(\Gamma_{-H}^R)$ and $\alpha_z \in H^{-1/2}(\Gamma_{-H}^R)$. However, to justify this second step, we adopt the following regularity assumption:

**Assumption 7.2.**

$$ \alpha_x, \alpha_z \in L^\infty(\Gamma_{-H}^R), \quad \exists m > 0, |\alpha_x(x)| \ge m \text{ for almost all } x \in \Gamma_{-H}^R. $$

In view of (54) and having in mind that $\theta^\delta$ is compactly supported in $(-R, R)$, a first simple method consists in computing directly $\theta^\delta$ as

$$ \theta^\delta(x) = \frac{1}{\alpha_x} \int_{-R}^{x} \alpha_z(t) \, dt. $$

A second method consists in using a weak formulation of (54). Let us denote $I = (-R, R)$ and consider the problem: find $\theta \in L^2(I)$ such that

$$ b(\theta, v) = \ell(v), \quad \forall v \in V_0(I), \qquad (58) $$
where the space $V_0(I)$ is defined by

$$V_0(I) = \{v \in H^1(I), v(-R) = 0\}$$

and where the sesquilinear form $b$ and the antilinear form $\ell$ are given by

$$b(\theta, v) := \int_I \alpha_x \theta \bar{v}' dx, \quad \ell(v) := - \int_I \alpha_z \bar{v} dx, \quad \forall (\theta, v) \in L^2(I) \times V_0(I). \quad (59)$$

We note that any solution to the weak formulation (58) satisfies the strong problem (54) in the sense of distributions. We now establish that the problem (58) is well-posed.

**Theorem 7.3.** Under Assumption 7.2, the problem (58) has a unique solution.

*Proof.* Let us apply the Brezzi-Nečas-Babuška Theorem (see for example [20]). Firstly, let us check the inf-sup condition, that is

$$\inf_{\theta \in L^2(I)} \sup_{v \in V_0(I)} \frac{|b(\theta, v)|}{\|\theta\|_{L^2(I)} \|v\|_{V_0(I)}} \ge \beta$$

for some $\beta > 0$. Indeed, we have

$$\sup_{v \in V_0(I)} \frac{|b(\theta, v)|}{\|v\|_{V_0(I)}} = \sup_{v \in V_0(I)} \frac{|\int_I \alpha_x \theta \bar{v}' dx|}{\|v\|_{V_0(I)}} = \|\alpha_x \theta\|_{L^2(I)} \ge m \|\theta\|_{L^2(I)},$$

where the second equality above is obtained by choosing $v$ such that

$$v(x) = \int_{-R}^{x} (\alpha_x \theta)(t) dt.$$

We hence obtain the result with $\beta = m$. Secondly, it is easily seen that if for some $v \in V_0(I)$ we have $b(\theta, v) = 0$ for all $\theta \in L^2(I)$, then $v = 0$. Indeed, if

$$\int_I \alpha_x \theta \bar{v}' dx = 0, \quad \forall \theta \in L^2(I),$$

then $\alpha_x v' = 0$, and given assumption 7.2 we get $v' = 0$, hence $v = 0$. $\square$

## 7.2 The acoustic model

Now let us consider the acoustic model. Again mimicking [3] or [9], we show that the first order approximation on $\Gamma_{-H}^R$ is $\partial_z u = 0$ while the second order approximation is

$$\partial_z u = \partial_x(\theta \partial_x u) + (\omega^2/c^2)\theta u \quad \text{on} \quad \Gamma_{-H}^R, \quad (60)$$
where $\|\theta\|_{L^\infty(-R,R)}$ is supposed to be small compared to the wavelength of the acoustic waves, that is $\lambda = 2\pi c/\omega$. In view of (33) we have

$$ \left\{ \begin{array}{rcll} \Delta u + (\omega^2/c^2)u & = & 0 & \text{in } \Omega_0^R, \\ u & = & 0 & \text{on } \Gamma_0^R, \\ -\partial_x u & = & T_-^a u - 2\partial_x \mathcal{G}^a(\cdot, N) & \text{on } \Sigma_{-R}, \\ \partial_x u & = & T_+^a u & \text{on } \Sigma_R. \end{array} \right. \qquad (61) $$

and (52), where the operators $T_\pm^a$ are defined by (31) while $\mathcal{G}^a(\cdot, N)$ is given by (32). We consider here the Hilbert spaces $V$, $M$ and $O$ as well as the operators $A$, $B$ and $C$ specified in section 6.2. We also define $L_1 \in M$ satisfying (53) and $G^a \in M$ satisfying (56) for $\mathcal{G}^a$ replaced by $\mathcal{G}^a$, $F_1 = L_1 + G^a$ and $F_2 = 0$, so that $F = (L_1 + G^a, 0) \in M \times O$. We easily check the following result.

**Lemma 7.4.** The function $u \in H^1(\Omega_0^R)$ satisfies the problem (61) (52) if and only if $u$ satisfies $Au = F$.

The above Lemma invites us again to directly apply all the results of section 5. In particular, the algorithm 5.2 enables us to obtain the Morozov solution $u^\delta \in H^1(\Omega_0^R)$ in the presence of some noisy data $F_1^\delta = L_1^\delta + G^a \in M$ such that (37) is satisfied.

In order to compute an approximated $\theta^\delta$ function in the presence of noise starting from the function $u^\delta$ obtained in the first step and reusing the definitions of $\alpha$, $\alpha_x$ and $\alpha_z$ given by (57), we wish to retrieve $\theta^\delta$ from the equation (60) in a second step. The weak formulation (58) is modified as follows: find $\theta \in L^2(I)$ such that

$$ b(\theta, v) + c(\theta, v) = \ell(v), \quad \forall v \in V_0(I), \qquad (62) $$

where the sesquilinear form $b$ and the antilinear form $\ell$ are defined by (59), while the sesquilinear form $c$ is given by

$$ c(\theta, v) := -(\omega^2/c^2) \int_I \alpha \theta \bar{v} \, dx, \quad \forall (\theta, v) \in L^2(I) \times V_0(I). $$

We now establish that the problem (62) is well-posed in the Fredholm sense.

**Theorem 7.5.** Under Assumption 7.2, if the problem (62) has at most one solution, such solution exists.

*Proof.* With the help of the Riesz theorem, we define the operators $B, C : L^2(I) \to V_0(I)$ such that

$$ (T\theta, v)_{V_0(I)} = t(\theta, v), \quad \forall (\theta, v) \in L^2(I) \times V_0(I), $$

where $T$ is either operator $B, C$ and $t$ is the corresponding sesquilinear form $b, c$. Since the sesquilinear form $b$ satisfies the assumptions of the Brezzi-Nečas-Babuška Theorem (see the proof of Theorem 7.3), the operator $B$ is
an isomorphism. Let us check that the operator $C$ is compact. To proceed, we set $u_{\theta} := C\theta \in V_0(I)$ for some $\theta \in L^2(I)$. Since from Assumption 7.2 we have $\alpha_x \in L^\infty(I) \subset L^2(I)$, we infer $\alpha \in H^1(I) \subset L^\infty(I)$. Hence $f := -(\omega^2/c^2)\alpha\theta \in L^2(I)$. We observe that the function $u_{\theta}$ is the unique solution in $V_0(I)$ to the weak formulation:

$$ \int_I u'_{\theta} \bar{v}' dx = \int_I f \bar{v} dx, \quad \forall v \in V_0(I). $$

A straightforward regularity result shows that $\theta \mapsto u_{\theta}$ is continuous from $L^2(I)$ to $H^2(I) \cap V_0(I)$, which implies that $C$ is compact. The conclusion results from the Fredholm alternative. $\square$

## 8 About the noisy data in the acoustic case

As it can be seen in the Algorithm 5.2 to compute the Morozov solution $u^\delta$, the amplitude of noise $\delta$ of the data $F$ is required in all cases. In the gravity case, such amplitude of noise coincides with the amplitude of noise on $F_2$, which happens to be, in view of (45), the amplitude of noise on the free surface perturbation $\eta$ multiplied by $(g/\omega)$. In the acoustic case, $\delta$ coincides with the amplitude of noise affecting $F_1$, that is the amplitude of noise on $L_1$. Estimating such amplitude of noise is more complicated than in the gravity case. Indeed, the volumic function $L_1$ is defined from the Neumann data $f_1$ by the weak formulation (53). In view of (52), the amplitude of noise on $f_1$ is given by multiplying the amplitude of noise on the free surface perturbation $\eta$ by $\omega$. But “converting” the amplitude of noise on $f_1$ into the amplitude on $L_1$ is not as simple. Let us denote

$$ \left\{ \begin{array}{ll} S: L^2(\Gamma_0^R) & \rightarrow M = \{\mu \in H^1(\Omega_0^R), \mu|_{\Gamma_{-H}^R} = 0\} \\ f & \mapsto L, \end{array} \right. \qquad (63) $$

where $L$ is the unique solution in $M$ of the weak formulation:

$$ \int_{\Omega_0^R} \nabla L \cdot \nabla \bar{\mu} \,dx dz = \int_{\Gamma_0^R} f \bar{\mu} \,dx, \quad \forall \mu \in M, \qquad (64) $$

which is equivalent to the strong problem: find $L \in M$ such that

$$ \left\{ \begin{array}{ll} \Delta L = 0 & \text{in } \Omega_0^R, \\ \partial_z L = f & \text{on } \Gamma_0^R, \\ L = 0 & \text{on } \Gamma_{-H}^R, \\ \partial_x L = 0 & \text{on } \Sigma_{\pm R}. \end{array} \right. \qquad (65) $$

If $d$ denotes the amplitude of noise on $f$, a simple way of estimating the resulting amplitude of noise $\delta$ on $L$ is

$$ \delta = \|S\|_d, \quad \|S\| := \sup_{f \in L^2(\Gamma_0^R), f \neq 0} \frac{\|Sf\|_M}{\|f\|_{L^2(\Gamma_0^R)}}. \qquad (66) $$
It happens that the operator norm $\|S\|$ is easily determined.

**Proposition 5.** We have $\|S\| = \sqrt{H}$ when $\theta = 0$.

*Proof.* Let us introduce the operator $T: L^2(\Gamma_0^R) \to L^2(\Gamma_0^R)$ defined by $T = \gamma \circ S$, where $\gamma: M \to L^2(\Gamma_0^R)$ is the trace operator. Let us check that the operator $T$ is self-adjoint, compact, positive and injective. Indeed, for $f, g \in L^2(\Gamma_0^R)$ we have

$$ (Tf, g)_{L^2(\Gamma_0^R)} = \int_{\Gamma_0^R} (Sf) \partial_z(\overline{Sg}) dx = \int_{\partial\Omega_0^R} (Sf) \partial_\nu(\overline{Sg}) ds \\ = \int_{\Omega_0^R} \nabla(Sf) \cdot \nabla(\overline{Sg}) dxdz = (f, Tg)_{L^2(\Gamma_0^R)}, \quad (67) $$

which proves that $T$ is self-adjoint. That $S$ is bounded and $\gamma$ is compact implies that the operator $T = \gamma \circ S$ is compact. Positivity and injectivity are straightforward consequences of (67). The spectral theorem implies then that the spectrum of $T$ is formed by a sequence of real and non increasing eigenvalues $\rho_n > 0$, $n \in \mathbb{N}$, such that $\rho_n \to 0$ when $n \to +\infty$, and there exists a complete orthonormal basis of $L^2(\Gamma_0^R)$ formed with the help of the corresponding eigenfunctions $e_n$. This implies that

$$ \|S\|^2 = \sup_{f \in L^2(\Gamma_0^R), f \neq 0} \frac{\|Sf\|_M^2}{\|f\|_{L^2(\Gamma_0^R)}^2} = \sup_{f \in L^2(\Gamma_0^R), f \neq 0} \frac{(Tf, f)_{L^2(\Gamma_0^R)}}{\|f\|_{L^2(\Gamma_0^R)}^2} = \rho_0, $$

since $\rho_0$ is the largest eigenvalue of $T$. An easy computation shows that

$$ \rho_0 = H, \quad \rho_n = \frac{2R}{n\pi} \tanh\left(\frac{n\pi H}{2R}\right), \quad n \ge 1, $$

and for $x \in (-R, R)$, we can choose

$$ e_0(x) = \frac{1}{\sqrt{2R}}, \quad e_n(x) = \frac{1}{\sqrt{R}} \cos\left(\frac{n\pi}{2R}(x+R)\right), \quad n \ge 1. $$

We have in particular $\rho_0 = H$, which completes the proof. $\square$

<figure><img src="image_2.png" /></figure>

The estimate (66) giving the amplitude $\delta$ of the noise on the data $L$ in $M$ from the amplitude of noise on the data $f$ in $L^2(\Gamma_0^R)$, in view of the definition of $\|S\|$, is likely to overestimate $\delta$. This is why we propose an alternative way of computing $\delta$ based on a probabilistic model of the measurements $f$. Let us assume that $f$ is given on $\Gamma_0^R$ by

$$ f = \sum_{k=1}^{N} f_k \chi_k^N. $$

For $k = 1, \dots, N$, the functions $\chi_k^N$ are indicator functions of intervals having the same length $2R/N$, which do not overlap and the union of which coincides with $\Gamma_0^R$. The complex random variables $f_k$, $k = 1, \dots, N$, satisfy the following discrete Gaussian white noise assumptions:
**Assumption 8.1.** The $f_k$, $k = 1, \dots, N$, are independent and identically distributed complex random variables such that their real and imaginary parts follow the normal distribution $N(0, \sigma^2/2)$, where $\sigma > 0$, which implies that

$$ \mathbb{E}[f_k] = 0, \quad \mathbb{E}[f_k \bar{f}_l] = \sigma^2 \delta_{kl}, \quad k, l = 1, \dots, N. $$

Instead of defining $\delta$ by (66), we can set

$$ \delta = r_N d, \quad r_N = \sqrt{\mathbb{E}\left[\frac{\|Sf\|_M^2}{\|f\|_{L^2(\Gamma_0^R)}^2}\right]}. \quad (68) $$

Such definition of $\delta$ seems more realistic: it accounts for the fact that the measurements are discrete in practice and that the error between the noisy and the exact data could reasonably satisfy Assumption 8.1. It is clear that $r_N \le \|S\|$, which implies that the resulting amplitude $\delta$ is smaller in the probabilistic approach than in the deterministic one. More precisely, the number $r_N$ is given by

**Theorem 8.2.**

$$ r_N = \sqrt{\frac{1}{2R} \sum_{k=1}^{N} \|S\chi_k^N\|_M^2}, $$

where we recall that $\| \cdot \|_M$ is the semi-norm in $H^1(\Omega_0^R)$.

The proof of Theorem 8.2 is omitted since it is very similar to that of Theorem 9.3 given in the appendix for a slightly simpler geometry. Besides, for this simpler geometry, in the appendix we also prove that $r_N \to 0$ when $N \to +\infty$. Theorem 8.2 offers a practical way of computing $r_N$: it requires to solve the weak formulation (64) $N$ times for $f = \chi_k^N$, $k = 1, \dots, N$.

# 9 Numerical experiments

## 9.1 The tsunami identification

All the numerical computations of the article are performed with the help of XLiFE++ Library [27]. We choose $c = 20$, $g = 10$, $H = 1$, which corresponds to $\omega_g \sim 3.16$ and $\omega_a = 20$. In this section, $R = 4$.

### 9.1.1 The gravity model

Let us begin with the tsunami identification using the gravity model, having in mind that the data are generated by the complete model. Starting from a function $\chi$ on $\Gamma_{-H}^R$ which represents the tsunami and that will be the reference function we try to retrieve, we compute the free surface perturbation
$\eta$ on $\Gamma_0^R$ by solving the complete forward problem (13) in the domain $\Omega_0^R$ with the help of a Finite Element Method. The mesh size corresponds to a subdivision of the segment $\Gamma_0^R$ into intervals of size $h = 0.01$, the finite elements consisting of classical P2 triangles. The infinite sums which define the Dirichlet-To-Neumann operators (12) are truncated to a finite number of terms which include all the propagating modes and at least 5 evanescent ones. Once the free surface perturbation $\eta$ has been computed, we artificially compute some noisy data $\eta^\delta$ on $\Gamma_0^R$ such that

$$ \eta^{\delta} = \eta + \rho b, $$

where $b$ is a function defined from its values at each discretization point of $\Gamma_0^R$. Such values are random complex numbers, the real and imaginary parts of which follow a uniform distribution $U(0, 1)$. The real $\rho > 0$ is calibrated such that

$$ \frac{\|\eta^{\delta} - \eta\|_{L^2(\Gamma_0^R)}}{\|\eta\|_{L^2(\Gamma_0^R)}} = \sigma, $$

where $\sigma$ is a prescribed relative noise amplitude. Multiplying the absolute amplitude of noise $\|\eta^{\delta} - \eta\|_{L^2(\Gamma_0^R)}$ by $g/\omega$ provides $\delta$. The artificial noisy data $\eta^\delta$ constitutes the input data of the Tikhonov-Morozov inverse method described in section 6.1, which lastly enables us to compute the retrieved tsunami function $\chi^\delta$. In the absence of noise, the retrieved tsunami function is given by $\chi_\varepsilon = (i/\omega)\partial_z u_\varepsilon$, where $u_\varepsilon$ is the solution to problem (48) for a very small $\varepsilon$. In order to discretize the problems (50) for all $n \in \mathbb{N}$, we again use a P2 Finite Element Method to compute the triple $(u_n, \lambda_n, h_n)$ and a mesh having the same size as for the forward problem. In the Figure 1, we compare the retrieved tsunami function $\chi^\delta$ and the actual tsunami function $\chi$ on $\Gamma_0^R$, for $\omega = 3$ and various relative amplitudes of noise $\sigma$ on $\eta$, namely $\sigma = 0$ (exact data), $\sigma = 1\%$, $\sigma = 5\%$ and $\sigma = 10\%$. Note that the value of $\omega$ is chosen close to $\omega_g$, that is the typical frequency of gravity waves, which explains the good reconstruction results for small $\sigma$.

### 9.1.2 The acoustic model

We now consider the tsunami identification with the help of the acoustic model. In the Figure 2, the retrieved function $\chi^\delta$ is compared to the exact function $\chi$, for $\omega = 20$ and the different relative amplitudes of noise on $\eta$, namely $\sigma = 0$ (exact data), $\sigma = 1\%$, $\sigma = 5\%$ and $\sigma = 10\%$. The value of $\omega$ is close to $\omega_a$, that is the typical frequency of acoustic waves. We use the same Finite Element Method and the same mesh size as in the gravity case. In the acoustic model, the amplitude of noise $\delta$ is obtained by multiplying $\|\eta^\delta - \eta\|_{L^2(\Gamma_0^R)}$ by $\omega$ and $r_N$, where $r_N$ is given by (68) with $N = 100$.
<figure><img src="image_3.png" /><figcaption>Figure 1: Gravity case for $\omega = 3$. Comparison between $\chi^\delta$ (dashed line) and $\chi$ (continuous line) for various amplitudes of noise. Top left: $\sigma = 0$. Top right: $\sigma = 1\%$. Bottom left: $\sigma = 5\%$. Bottom right: $\sigma = 10\%$.</figcaption></figure>

**Remark 4.** In the Figures 1 and 2, we notice that the error between the retrieved tsunami $\chi^\delta$ and the exact one $\chi$ is quite large at the endpoints of $\Gamma_{-H}^R$, which correspond to the transition points between the part of the boundary $\partial\Omega^R$ on which we have data and the part $\Gamma_{-H}^R$ on which we have no data. As shown by the analysis conducted in [8], we explain this phenomenon by the fact that the smoothness of the solution to the mixed formulation of the Tikhonov regularization is poorer at those transition points. Even if such error decreases when the mesh size is refined, we have decided to apply an additional regularization to $\chi^\delta$ on $\Gamma_{-H}^R$. From the a priori assumption that $\chi$ is compactly supported in $I = (-R, R)$, we define $\tilde{\chi}^\delta$ as the solution of the one dimensional weak formulation: find $\tilde{\chi}^\delta \in H_0^1(I)$ such that

$$ \int_I \tilde{\chi}^\delta \bar{v} dx + \tilde{\varepsilon} \int_I (\tilde{\chi}^\delta)' \bar{v}' dx = \int_I \chi^\delta \bar{v} dx, \quad \forall v \in H_0^1(I). $$

Here, the regularization parameter is given by $\tilde{\varepsilon} = 0.001/k^2$, where $k = \omega/\sqrt{gH}$ in the gravity case and $k = \omega/c$ in the acoustic case. Such regularization was already applied in Figures 1 and 2, otherwise the errors at the endpoints would have been worse. In what follows the regularized function
<figure><img src="image_4.png" /><figcaption>Figure 2: Acoustic case for $\omega = 20$. Comparison between $\chi^\delta$ (dashed line) and $\chi$ (continuous line) for various amplitudes of noise. Top left: $\sigma = 0$. Top right: $\sigma = 1\%$. Bottom left: $\sigma = 5\%$. Bottom right: $\sigma = 10\%$.</figcaption></figure>

$\tilde{\chi}^\delta$ will be simply denoted $\chi^\delta$.

### 9.1.3 On the relevance of the Morozov choice

We wish now to test the relevance of the Morozov strategy to choose the regularization parameter $\varepsilon$ as a function of the amplitude of noise $\delta$. In the acoustic case, we have plotted in the Figure 3 the retrieved tsunami $\chi_\varepsilon^\delta = (i/\omega)\partial_z u_\varepsilon^\delta$ and the actual one $\chi$, where $u_\varepsilon^\delta$ is the solution to the problem (35) corresponding to the acoustic case, by using a value of $\varepsilon$ in (35) which seems relevant but differs from the Morozov choice, for two relative amplitudes of noise on $\eta$. Comparing the two bottom pictures of Figure 2 and the pictures of Figure 3, we realize that the Morozov choice for $\varepsilon$ is not necessarily the best in terms of the error between the retrieved and the exact solution. To be more precise, in Figure 4 we have plotted the error between the retrieved tsunami $\chi_\varepsilon^\delta$ and the actual one $\chi$ as a function of $\varepsilon$. In Figure 4 we also compare the deterministic and probabilistic choices of $\delta$ described in section 8, which are given by (66) and (68), respectively. Note that the $L^2$ error is limited to the subpart $\Gamma_{-H}^{R'}$ of $\Gamma_{-H}^R$ with $R' = 0.9R$, since the reconstruction
<figure><img src="image_5.png" /><figcaption>Figure 3: Acoustic case for $\omega = 20$. Comparison between $\chi_{\varepsilon}^{\delta}$ (dashed line) and $\chi$ (continuous line) for two amplitudes of noise and well-tuned choices for $\varepsilon$ (see Figure 4). Left: $\sigma = 5\%$. Right: $\sigma = 10\%$.</figcaption></figure>

close to the endpoints of $\Gamma_{-H}^{R}$ is bad and would prevent us from presenting a clear analysis (see Remark 4). Classically, we observe that in the absence

<figure><img src="image_6.png" /><figcaption>Figure 4: Error $\|\chi_{\varepsilon}^{\delta} - \chi\|_{L^2(\Gamma_{0}^{R'})}$ between the solution of problem (35) in the acoustic case and the exact solution, as a function of $\varepsilon$. The blue curve corresponds to uncontaminated data $\sigma = 0\%$ while the red one corresponds to noisy data with $\sigma = 5\%$. The vertical line on the right corresponds to the Morozov choice for $\delta$ following (66) while the other vertical dashed lines correspond to the Morozov choice for $\delta$ following (68) with $N = 100$, $N = 400$ and $N = 800$, from the right to the left.</figcaption></figure>

of noise ($\delta = 0$) the error $\|\chi_{\varepsilon} - \chi\|_{L^2(\Gamma_{0}^{R'})}$ increases with $\varepsilon$, except for very small values of $\varepsilon$, due to the finite element discretization. In the presence of noise, the error $\|\chi_{\varepsilon}^{\delta} - \chi\|_{L^2(\Gamma_{0}^{R'})}$ first decreases with $\varepsilon$ and then increases with $\varepsilon$, which justifies to look for an optimal value of $\varepsilon$. We observe that the value of $\varepsilon$ given by the Morozov choice is better by using the probabilistic strategy (68) than the deterministic strategy (66). More precisely the best
value of $N$ is 400, which roughly corresponds to the number of discretization points of $\Gamma_0^R$, that is the number of artificial measurements, which is expected.

### 9.1.4 Influence of the chosen model

In order to compare the quality of the identification produced by the different models (gravity, acoustics and complete), in Figure 5 we have plotted the error between the retrieved tsunami function $\chi^\delta$ using the Tikhonov-Morozov strategy and the actual one $\chi$ as a function of the frequency $\omega$ on the subpart $\Gamma_{-H}^{R'}$. When data are uncontaminated, the best model for the identification

<figure><img src="image_7.png" /><figcaption>Figure 5: Error $\|\chi^\delta - \chi\|_{L^2(\Gamma_{-H}^{R'})}$ as a function of $\omega$ using different models: gravity (green), acoustic (red) and complete (blue). Left: uncontaminated data ($\sigma = 0\%$). Right: noisy data ($\sigma = 5\%$)</figcaption></figure>

seems to be the complete one for all frequencies, which is expected since the data were produced by the complete model. The results obtained with the gravity model are good at a very low frequency but strongly deteriorates when $\omega$ increases, which was expected in view of Remark 2. In the presence of noisy data, the best model is again the complete one except at a very low frequency, for which the acoustic one is better. Note that the gravity model is competitive at frequencies around 10, strongly deteriorates at high frequencies, and more surprinsigly also at very low frequencies.

### 9.1.5 The case of a variable sea bottom

We complete this section by showing that our strategy to identify a tsunami can easily be extended to a non-flat sea bottom, contrary to the other methods mentioned in the introduction, which rely on integral transforms with respect to the horizontal variable. It suffices to use two different left and right Dirichlet-To-Neumann operators, those operators being defined with the help of the guided modes related to the left and right half-waveguides, respectively. This is illustrated by Figure 6 in the acoustic case, with $\omega = 20$ and $\sigma = 1\%$.
<figure><img src="image_8.png" /><figcaption>Figure 6: Non flat bottom, acoustic case for $\omega = 20$. Left: Geometry of the sea bottom. Right: Comparison between $\chi^\delta$ (dashed line) and $\chi$ (continuous line) for $\sigma = 1\%$.</figcaption></figure>

## 9.2 Bathymetry

In this section, $R = 1$ and $H = 0.5$, the other parameters being unchanged. Let us begin with the bathymetry problem using the gravity model, the data still being obtained with the complete model. Starting from a reference function $\theta$ which governs the profile $\Gamma_{-H+\theta}^R$ of the sea bottom, we compute the free surface perturbation $\eta$ on $\Gamma_0^R$ by solving the forward problem (19) in the domain $\Omega_\theta^R$ with the help of the same Finite Element Method as for the tsunami identification and the same mesh size. The source point $N$ is chosen as $(-1.1, -0.45)$. Then we artificially perturb $\eta$ to generate a noisy data $\eta^\delta$ on $\Gamma_0^R$ by using the same procedure as for the tsunami identification. Lastly we use the two steps strategy based on the approximate second-order boundary condition on $\Gamma_{-H}^R$ and described in section 7.1. It is important to note that the data are obtained by solving a forward problem based on the true non-flat geometry $\Gamma_{-H+\theta}^R$ while the inverse problem is solved in a flat geometry $\Gamma_{-H}^R$. We finally obtain a function $\theta^\delta$ on $\Gamma_{-H}^R$ in the presence of noisy data. Note that in the second step, such $\theta^\delta$ is obtained with the help of a finite element discretization of the one dimensional weak formulation (58). In order to discretize the spaces $L^2(I)$ and $V_0(I)$, we use the same finite dimensional space based on $P1$ finite elements. In Figure 7 we compare the retrieved function $\theta^\delta$ and the reference one $\theta$, both with uncontaminated data ($\sigma = 0$) and with noisy data ($\sigma = 1\%$), for a frequency $\omega = 3$ compatible with the gravity model. Let us secondly consider the bathymetry problem using the acoustic model. The two step strategy is now described by section 7.2. In the second step, the computation of the function $\theta^\delta$ consists in a finite element discretization of the weak formulation (62) which is the same as in the gravity case. In Figure 8 we compare the retrieved function $\theta^\delta$ and the true one $\theta$, both with uncontaminated data ($\sigma = 0$) and with noisy data ($\sigma = 1\%$), for frequency $\omega = 25$. Unsurprisingly, the identification results
<figure><img src="image_9.png" /><figcaption>Figure 7: Gravity case for $\omega = 3$. Comparison between $\theta^\delta$ (dashed line) and $\theta$ (continuous line). Left: uncontaminated data ($\sigma = 0$). Right: noisy data ($\sigma = 1\%$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8: Acoustic case for $\omega = 25$. Comparison between $\theta^\delta$ (dashed line) and $\theta$ (continuous line). Left: uncontaminated data ($\sigma = 0$). Right: noisy data ($\sigma = 1\%$).</figcaption></figure>

for bathymetry are more sensitive with respect to $\sigma$ than those obtained for the tsunami identification, which is due to the fact that the former is a geometric nonlinear inverse problem while the latter is linear. In order to compare the three different models with respect to the frequency $\omega$, in Figure 9 we have plotted the error between the retrieved function $\theta^\delta$ and the actual one as a function of the frequency $\omega$. We observe that the results are better for the complete and acoustic models than for the gravity one whatever the frequency is, which can be explained by Remark 2. In the gravity case, and as for the tsunami identification problem, the results are satisfactory for small $\omega$ but strongly deteriorate when $\omega$ increases. This phenomenon also occurs for the complete/acoustic model, for $\omega$ larger than 30: this is probably due to the fact that the approximation of the true boundary condition at the bottom by a second-order impedance boundary condition becomes false
<figure><img src="image_11.png" /><figcaption>Figure 9: Error $||\theta^\delta - \theta||_{L^2(\Gamma_{-H}^R)}$ as a function of $\omega$ using different models: gravity (green), acoustic (red) and complete (blue). Left: uncontaminated data ($\sigma = 0\%$). Right: noisy data ($\sigma = 1\%$).</figcaption></figure>

when the wavelength $\lambda := 2\pi c/\omega$ is not large compared to the amplitude of $\theta$. It is important to note that in the numerical results related to bathymetry, the assumption 7.2 seems to be satisfied.

From the numerical results obtained for the tsunami problem and for the bathymetry problem, we can conclude that these two inverse problems are approximately solvable using the gravity model only at a very low frequency (typically $\omega \simeq \sqrt{g/H}$) while they are approximately solvable at any (not too large) frequency for the acoustic/complete models. Note that the identification results were slightly better using the complete model than the acoustic one, but the reader has to keep in mind that the artificial data were obtained by using the complete model in any case.

## Appendix : a Laplace type problem with random Neumann boundary conditions

In order to study the limit of $r_N$ defined by (68) when $N \to +\infty$, instead of considering the problem (65) we move to a slightly simpler problem in the two dimensional disc $D$ centered at 0 and of radius $R$: for some boundary data $f$, find $L$ such that

$$ \left\{ \begin{array}{ll} -\Delta L + L = 0 & \text{in } D, \\ \partial_\nu L = f & \text{on } \partial D. \end{array} \right. \qquad (69) $$

The problem (69) is obviously well-posed in $H^1(D)$ for $f \in L^2(\partial D)$. We can hence define the operator

$$ \left\{ \begin{array}{ll} S : L^2(\partial D) \to H^1(D) \\ f \mapsto L, \end{array} \right. \qquad (70) $$
such that $L$ is the solution to problem (69), as well as the operator $T = \gamma \circ S : L^2(\partial D) \to L^2(\partial D)$, where $\gamma : H^1(D) \to L^2(\partial D)$ is the trace operator. A similar proposition as Proposition 5 is now obtained.

**Proposition 6.** We have

$$ \|S\| = \sqrt{\frac{I_0(R)}{I'_0(R)}}, $$

where for $n \in \mathbb{N}$, the $I_n(r)$ are the bounded solutions (up to a multiplicative constant) to the modified Bessel equations

$$ r^2 \frac{d^2 u_n}{dr^2} + r \frac{du_n}{dr} - (n^2 + r^2)u_n = 0, \quad r \in (0, R). $$

**Remark 5.** A comprehensive analysis of the functions $I_n(r)$ may be found in [1].

*Proof of Proposition 6.* By using the same arguments as in the proof of Proposition 5, we show that the operator $T$ is self-adjoint, compact, positive and injective. The spectrum of $T$ is formed by a sequence of nonincreasing eigenvalues $\rho_n > 0$, $n \in \mathbb{N}$, such that $\rho_n \to 0$ when $n \to +\infty$, the corresponding eigenfunctions $e_n$ constituting a complete orthonormal basis of $L^2(\partial D)$. As in the proof of Proposition 5, we also have $\|S\|^2 = \rho_0$. An easy computation shows that

$$ \rho_n = \frac{I_n(R)}{I'_n(R)}, \quad e_n(\theta) = \frac{1}{\sqrt{2\pi R}} e^{in\theta}, \quad n \in \mathbb{N}, $$

which implies the stated result.

<figure><img src="image_12.png" /></figure>

**Remark 6.** We observe that the operator $T$ is a Hilbert-Schmidt operator (see for example [13] for a definition), since

$$ \sum_{n=0}^{+\infty} \|Te_n\|_{L^2(\partial D)}^2 = \sum_{n=0}^{+\infty} \rho_n^2 < +\infty. $$

Indeed, by using the asymptotics of the modified Bessel function given by [1], we have

$$ \rho_n \underset{n \to +\infty}{\sim} \frac{R}{n}, $$

which enables us to conclude.

In the following theorem, we prove that the problem (69) is also well-posed in $L^2(D)$ for $f \in H^{-3/2}(\partial D)$.

**Theorem 9.1.** For $f \in H^{-3/2}(\partial D)$, the problem (69) has a unique solution in $L^2(D)$.
Proof. Let us introduce the space $V := \{v \in H^2(D), \partial_\nu v|_{\partial D} = 0\}$ and the operator

$$ \left\{ \begin{array}{l} A: V \rightarrow L^2(D) \\ v \mapsto -\Delta v + v. \end{array} \right. $$

From a standard regularity result for elliptic equations (see for example [13]), the continuous operator $A$ is an isomorphism. This implies that the adjoint operator $A^* : L^2(D) \rightarrow V^*$, where $V^*$ is the dual space of $V$, is also an isomorphism (see again [13]). This means that for all $L^* \in V^*$, there exists a unique $u \in L^2(D)$ such that $A^*u = L^*$. Such $u$ satisfies

$$ \langle A^*u, v \rangle_{V^*, V} = (u, Av)_{L^2(D)} = \langle L^*, v \rangle_{V^*, V}, \quad \forall v \in V. $$

If we choose $L^* \in V^*$ such that

$$ \langle L^*, v \rangle_{V^*, V} = \langle f, v \rangle_{H^{-3/2}(\partial D), H^{3/2}(\partial D)}, \quad \forall v \in V, $$

we conclude that for $f \in H^{-3/2}(\partial D)$, there exists a unique $u \in L^2(D)$ such that

$$ \int_D u(-\Delta v + v) dx = \langle f, v \rangle_{H^{-3/2}(\partial D), H^{3/2}(\partial D)}, \quad \forall v \in V. \quad (71) $$

Let us check that such function $u$ satisfies problem (69). Taking $v = \varphi \in C_0^\infty(D)$, we first obtain that $-\Delta u + u = 0$ in $D'(D)$. Hence $u, \Delta u \in L^2(D)$, so that for all $v \in V$, we have the Green formula (see [31]):

$$ \int_D u(-\Delta v + v) dx = \int_D (-\Delta u + u)v dx - \langle u, \partial_\nu v \rangle_{H^{-1/2}(\partial D), H^{1/2}(\partial D)} + \langle \partial_\nu u, v \rangle_{H^{-3/2}(\partial D), H^{3/2}(\partial D)}. $$

In view of (71), we obtain that

$$ \langle \partial_\nu u - f, v \rangle_{H^{-3/2}(\partial D), H^{3/2}(\partial D)} = 0, \quad \forall v \in V. $$

We conclude that $\partial_\nu u = f$ in $H^{-3/2}(\partial D)$. $\square$

For all $y \in \partial D$, since dimension is 2 we have $\delta_y \in H^{-3/2}(\partial D)$, so that for $f = \delta_y$, problem (69) has a unique solution in $L^2(D)$, which is denoted $G_y = G(\cdot, y)$.

The operator $T: L^2(\partial D) \rightarrow L^2(\partial D)$ can then be seen as an integral operator, defined for all $f \in L^2(\partial D)$, by

$$ (Tf)(x) = \int_{\partial D} G(x, y)f(y) ds, \quad \forall x \in \partial D. \quad (72) $$

**Remark 7.** From Remark 6, since $T$ is a Hilbert-Schmidt operator, we have that $G \in L^2(\partial D \times \partial D)$ (see for example [13]).
Like in section 8, let us assume that $f$ is given on $\partial D$ by

$$f = \sum_{k=1}^{N} f_k \chi_k^N. \qquad (73)$$

For $k = 1, \dots, N$, the functions $\chi_k^N$ are indicator functions of intervals having the same length $2\pi R/N$, which do not overlap and the union of which coincides with $\partial D$. The random variables $f_k$, $k = 1, \dots, N$, satisfy the following discrete Gaussian white noise assumptions:

**Assumption 9.2.** The $f_k$, $k = 1, \dots, N$, are independent and identically distributed real random variables following the normal distribution $N(0, \sigma^2)$, where $\sigma > 0$, which implies

$$\mathbb{E}[f_k] = 0, \quad \mathbb{E}[f_k f_l] = \sigma^2 \delta_{kl}, \quad k, l = 1, \dots, N.$$

We again introduce

$$r_N = \sqrt{\mathbb{E}\left[\frac{\|Sf\|_{H^1(D)}^2}{\|f\|_{L^2(\partial D)}^2}\right]}. \qquad (74)$$

We now derive a simple expression for $r_N$.

**Theorem 9.3.** The number $r_N$ is given by

$$r_N = \sqrt{\frac{N}{2\pi R}} \|S\chi_1^N\|_{H^1(D)}.$$

*Proof.* The weak formulation of problem (69) implies that

$$\|Sf\|_{H^1(D)}^2 = \int_{\partial D} f(Tf) ds, \quad \forall f \in L^2(\partial D).$$

In view of (72), we deduce the concrete expression

$$\|Sf\|_{H^1(D)}^2 = \int_{\partial D} \int_{\partial D} G(x, y) f(x) f(y) ds(x) ds(y).$$

The decomposition (73) yields

$$\|Sf\|_{H^1(D)}^2 = (G^N X, X)_{\mathbb{R}^N},$$

where the vector $X \in \mathbb{R}^N$ and the matrix $G^N \in \mathbb{R}^{N \times N}$ are defined by

$$X_k := f_k, \quad G_{kl}^N = \int_{\partial D} \int_{\partial D} G(x, y) \chi_k^N(x) \chi_l^N(y) ds(x) ds(y), \quad k, l = 1, \dots, N.$$
On the other hand,

$$ \|f\|_{L^2(\partial D)}^2 = \sum_{k,l=1}^{N} f_k f_l \int_{\partial D} \chi_k^N \chi_l^N ds. $$

From the definition of the $\chi_k^N$, we obtain

$$ \|f\|_{L^2(\partial D)}^2 = \frac{2\pi R}{N} \sum_{k=1}^{N} f_k^2 = \frac{2\pi R}{N} \|X\|_{\mathbb{R}^N}^2. $$

We conclude that

$$ r_N^2 = \frac{N}{2\pi R} \mathbb{E} \left[ \frac{(G^N X, X)_{\mathbb{R}^N}}{\|X\|_{\mathbb{R}^N}^2} \right]. $$

In view of Assumption 9.2 and Theorem 1 in [25], we have

$$ \mathbb{E} \left[ \frac{(G^N X, X)_{\mathbb{R}^N}}{\|X\|_{\mathbb{R}^N}^2} \right] = \frac{\mathbb{E} [(G^N X, X)_{\mathbb{R}^N}]}{\mathbb{E} [\|X\|_{\mathbb{R}^N}^2]} $$

We easily derive, again by using Assumption 9.2, that

$$ \mathbb{E} [\|X\|_{\mathbb{R}^N}^2] = N\sigma^2 $$

and

$$ \mathbb{E}[(G^N X, X)_{\mathbb{R}^N}] = \sum_{k,l=1}^{N} G_{kl}^N \mathbb{E}[f_k f_l] = \sigma^2 \sum_{k=1}^{N} G_{kk}^N = N\sigma^2 G_{11}^N. $$

In the last equality, we have used the symmetry of $\partial D$, which implies that $G_{kk}^N$ does not depend on $k = 1, \dots, N$. We conclude that

$$ r_N^2 = \frac{N}{2\pi R} G_{11}^N, $$

with

$$ G_{11}^N = \int_{\partial D} \int_{\partial D} G(x, y) \chi_1^N(x) \chi_1^N(y) ds(x)ds(y) = \|S\chi_1^N\|_{H^1(D)}^2, $$

which completes the proof. □

<figure><img src="image_13.png" /></figure>

An interesting question concerns the behaviour of $r_N$ when $N \to +\infty$. In order to study this limit, we need to specify the singularity of the kernel $G$ (Remark 7 is not sufficient to conclude).

**Lemma 9.4.** The function $G$ satisfies $G \in L^p(\partial D \times \partial D)$, for all $p \ge 1$.
Proof. Let us introduce the function $(x, y) \mapsto G^{\infty}(x, y) := -\ln|x - y|/\pi$ and the function $x \mapsto G_y^{\infty}(x) := G^{\infty}(x, y)$ for a given $y$. The latter satisfies $\Delta G_y^{\infty} = 0$ for $x \neq y$. We clearly have that $G^{\infty} \in L^p(\partial D \times \partial D)$ for all $p \ge 1$. Then it is enough to prove that $G - G^{\infty}$ is bounded in $\partial D \times \partial D$, which constitutes the remainder of the proof. Let us pick some $y \in \partial D$. From (71), the function $G_y \in L^2(D)$ satisfies:

$$ \int_D G_y(-\Delta v + v) dx = v(y), \quad \forall v \in V. \qquad (75) $$

On the other hand, by denoting $D_{\varepsilon} := D \setminus \overline{D(y, \varepsilon)}$, where $D(y, \varepsilon)$ is the disc centered at $y$ and of radius $\varepsilon > 0$,

$$ \int_D G_y^{\infty}(-\Delta v + v) dx = \lim_{\varepsilon \to 0} \int_{D_{\varepsilon}} G_y^{\infty}(-\Delta v + v) dx \\ = \lim_{\varepsilon \to 0} \left\{ \int_{D_{\varepsilon}} v(-\Delta G_y^{\infty} + G_y^{\infty}) dx - \int_{\partial D_{\varepsilon}} G_y^{\infty} \partial_{\nu} v ds + \int_{\partial D_{\varepsilon}} \partial_{\nu} G_y^{\infty} v ds \right\}, $$

where the unit normal $\nu$ is oriented outside $D_{\varepsilon}$. We have $\Delta G_y^{\infty} = 0$ in $D_{\varepsilon}$, and

$$ \lim_{\varepsilon \to 0} \int_{D_{\varepsilon}} v(-\Delta G_y^{\infty} + G_y^{\infty}) dx = \int_{D} G_y^{\infty} v dx. $$

We focus now on the limit when $\varepsilon \to 0$ of the two boundary terms. Since $\partial_{\nu} v = 0$ on $\partial D$, we have

$$ \int_{\partial D_{\varepsilon}} G_y^{\infty} \partial_{\nu} v ds = \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} G_y^{\infty} \partial_{\nu} v ds = -\frac{\ln \varepsilon}{\pi} \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} \partial_{\nu} v ds. $$

But

$$ \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} \partial_{\nu} v ds = \int_{\partial(D \cap D(y, \varepsilon))} \partial_{\nu} v ds = \int_{D \cap D(y, \varepsilon)} \Delta v dx. $$

By the Cauchy-Schwarz inequality, it follows that

$$ \left| \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} \partial_{\nu} v \, ds \right| \le \| \Delta v \|_{L^2(D \cap D(y, \varepsilon))} \sqrt{\int_{D \cap D(y, \varepsilon)} 1 \, dx}, $$

with

$$ \sqrt{\int_{D \cap D(y, \varepsilon)} 1 \, dx} \le \sqrt{\frac{\pi}{2}} \varepsilon. $$

Gathering the previous estimates, we conclude that there exists a constant $c > 0$ such that

$$ \left| \int_{\partial D_{\varepsilon}} G_y^{\infty} \partial_{\nu} v \, ds \right| \le c \varepsilon |\ln \varepsilon| \| \Delta v \|_{L^2(D)} \xrightarrow{\varepsilon \to 0} 0. $$
Concerning the second boundary term, we have

$$ \int_{\partial D_{\varepsilon}} \partial_{\nu} G_{y}^{\infty} v \, ds = \int_{\partial D_{\varepsilon} \cap \partial D} \partial_{\nu} G_{y}^{\infty} v \, ds + \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} \partial_{\nu} G_{y}^{\infty} v \, ds. $$

Considering the first term above, we observe that the function $g_y$ defined by

$$ g_y(x) := -\partial_{\nu} G_y^{\infty}(x) = \frac{1}{\pi} \frac{1}{|x-y|^2} (x-y) \cdot \nu_x, \quad \forall x \in \partial D \setminus \{y\}, $$

is continuous on $\partial D$. Indeed, let us assume that the boundary $\partial D$ is locally characterized by the equation $h(x) = 0$, where $h$ is an infinitely smooth function. For $x \in \partial D$ and $x$ is sufficiently close to $y$,

$$ h(y) = h(x) + \nabla h(x) \cdot (y-x) + \frac{1}{2} (y-x)^T \cdot \nabla^2 h(x) \cdot (y-x) + O(|x-y|^3), $$

where $\nabla h$ and $\nabla^2 h$ are the gradient and the Hessian of $h$, respectively, with $\nabla h(x) \neq 0$ for $x$ close to $y$. Since $h(x) = h(y) = 0$ and $\nu_x = \nabla h(x)/|\nabla h(x)|$, this yields

$$ g_y(x) = \frac{1}{2\pi} \frac{(x-y)^T}{|x-y|} \cdot \frac{\nabla^2 h(x)}{|\nabla h(x)|} \cdot \frac{(x-y)}{|x-y|} + O(|x-y|) \xrightarrow{x \to y} \frac{1}{2\pi} \tau_y^T \cdot \frac{\nabla^2 h(y)}{|\nabla h(y)|} \cdot \tau_y, $$

where $\tau_y$ is the tangential vector to $\partial D$ at point $y$. This proves that $g_y \in C^0(\partial D)$. Finally,

$$ \int_{\partial D_{\varepsilon} \cap \partial D} \partial_{\nu} G_y^{\infty} v \, ds \xrightarrow{\varepsilon \to 0} - \int_{\partial D} g_y v \, ds. $$

Considering the second term, we have

$$ g_y(x) = \frac{1}{\pi} \frac{1}{|x-y|^2} (x-y) \cdot \nu_x = -\frac{1}{\pi \varepsilon}, \quad \forall x \in \partial D_{\varepsilon} \cap \partial D(y, \varepsilon), $$

hence

$$ \int_{\partial D_{\varepsilon} \cap \partial D(y, \varepsilon)} \partial_{\nu} G_y^{\infty} v \, ds \xrightarrow{\varepsilon \to 0} v(y) \frac{1}{\pi \varepsilon} = v(y). $$

As a conclusion, the function $G_y^\infty$ satisfies the weak formulation

$$ \int_D G_y^\infty(-\Delta v + v) \, dx = \int_D G_y^\infty v \, dx - \int_{\partial D} g_y v \, ds + v(y), \quad \forall v \in V. \quad (76) $$

Subtracting the weak formulations (75) and (76), we obtain that the function $r_y := G_y - G_y^\infty \in L^2(D)$ satisfies the weak formulation:

$$ \int_D r_y(-\Delta v + v) \, dx = - \int_D G_y^\infty v \, dx + \int_{\partial D} g_y v \, ds, \quad \forall v \in V, $$
then $r_y$ satisfies the strong problem:

$$ \begin{cases} -\Delta r_y + r_y = -G_y^\infty & \text{in } D, \\ \partial_\nu r_y = g_y & \text{on } \partial D. \end{cases} $$

Since $G_y^\infty \in L^2(D)$ and $g_y \in L^2(\partial D)$, by a standard regularity result we get that $r_y \in H^{3/2}(D)$, then $r_y \in C^0(\bar{D})$, in particular $r_y|_{\partial D} \in C^0(\partial D)$. Hence $r_y$ is bounded on $\partial D$ for all $y \in \partial D$ and the bound $\|r_y\|_{L^\infty(\partial D)}$ does not depend on $y$ by rotation invariance. Finally, $G - G_y^\infty$ is bounded on $\partial D \times \partial D$, which completes the proof. $\square$

Now let us specify the behaviour of $r_N$ when the number $N$ of measurements becomes large.

**Theorem 9.5.** We have

$$ \lim_{N \to +\infty} N^{1/2-\epsilon} r_N = 0, \quad \forall \epsilon > 0. $$

*Proof.* We recall that

$$ r_N^2 = \frac{N}{2\pi R} \int_{\partial D} \int_{\partial D} G(x, y) \chi_1^N(x) \chi_1^N(y) ds(x)ds(y). $$

From Lemma 9.4, we have in particular that $G \in L^q(\partial D \times \partial D)$, for all $q > 2$. Using the Hölder inequality, with $p$ defined by $1/p + 1/q = 1$, we have

$$ \begin{aligned} r_N^2 &\le \frac{N}{2\pi R} \left( \int_{\partial D \times \partial D} |G(x, y)|^q ds(x)ds(y) \right)^{1/q} \\ &\quad \times \left( \int_{\partial D \times \partial D} |\chi_1^N(x)\chi_1^N(y)|^p ds(x)ds(y) \right)^{1/p} \\ &= \frac{N}{2\pi R} \|G\|_{L^q(\partial D \times \partial D)} \left( \frac{2\pi R}{N} \times \frac{2\pi R}{N} \right)^{1/p} \\ &= \|G\|_{L^q(\partial D \times \partial D)} \left( \frac{2\pi R}{N} \right)^{1-2/q}. \end{aligned} $$

That $q > 2$ enables us to conclude. $\square$

In Figure 10 below, we numerically illustrate both Theorem 9.3 and Theorem 9.5. More precisely, in the left part of Figure 10 we compare the value of $r_N$ defined by (74), which is computed with the help of Theorem 9.3, and its Monte-Carlo approximation given by the empirical mean

$$ (r_N^M)^2 = \frac{1}{M} \sum_{m=1}^{M} \frac{\|Sf^m\|_{H^1(D)}^2}{\|f^m\|_{L^2(\partial D)}^2}, \quad (77) $$

where the $f^m$, $m = 1, \dots, M$, form a sample of $M$ generated functions $f$ which satisfy Assumption 9.2. We can check the convergence of $r_N^M$ to $r_N$
for large M. In the right part of Figure 10, we have plotted the values of $r_N$ with respect to $N$ in a log-log scale, which exhibits a negative slope of 0.39 approximately. Such slope is quite close to the optimal value given by Theorem 9.5. Besides, the value of $\|S\|$, which majorates $r_N$ for all $N$ according to (74) and is given by Proposition 6, is such that $\log \|S\| = 0.175$. That $r_N \to 0$ when $N \to +\infty$ means in some sense that noise is vanishing

<figure><img src="image_14.png" /><figcaption>Figure 10: Left: Convergence of the Monte-Carlo estimates $r_N^M$ to $r_N$ when $M \to +\infty$. Right: Values of $\log r_N$ with respect to $\log N$ ($\log \|S\| = 0.175$).</figcaption></figure>

when the number of measurements tends to infinity. We note however that in practice $N$ is the number of measurements, which is finite and known.

## Acknowledgement

The authors are indebted to Lucas Chesnel for helpful discussions.

## References

[1] Milton Abramowitz and Irene A. Stegun. *Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables*. Applied mathematics series. Dover Publications, 1965.

[2] Valeria Bacchelli. Uniqueness for the determination of unknown boundary and impedance with the homogeneous robin condition. *Inverse Problems*, 25(1):015004, nov 2008.

[3] Aslanyürek Baslan, Haddar Houssem, and Şahintürk Hülya. Generalized impedance boundary conditions for thin dielectric coatings with variable thickness. *Wave Motion*, 48(7):681-700, 2011.

[4] Anne-Sophie Bonnet-Bendhia, Lynda Dahi, Eric Lunéville, and Vincent Pagneux. Acoustic diffraction by a plate in a uniform flow. *Mathematical Models and Methods in Applied Sciences*, 12(05):625-647, 2002.
[5] Laurent Bourgeois. A mixed formulation of quasi-reversibility to solve the Cauchy problem for Laplace's equation. *Inverse Problems*, 21(3):1087-1104, 2005.

[6] Laurent Bourgeois. Convergence rates for the quasi-reversibility method to solve the cauchy problem for laplace's equation. *Inverse problems*, 22(2):413, 2006.

[7] Laurent Bourgeois, Nicolas Chaulet, and Houssem Haddar. On simultaneous identification of the shape and generalized impedance boundary condition in obstacle scattering. *SIAM Journal on Scientific Computing*, 34(3):A1824–A1848, 2012.

[8] Laurent Bourgeois and Lucas Chesnel. On quasi-reversibility solutions to the cauchy problem for the laplace equation: regularity and error estimates. *ESAIM: M2AN*, 54(2):493-529, 2020.

[9] Laurent Bourgeois and Lucas Chesnel. Generalized impedance boundary conditions with vanishing or sign-changing impedance. working paper or preprint, September 2023.

[10] Laurent Bourgeois and Jérémi Dardé. A duality-based method of quasi-reversibility to solve the cauchy problem in the presence of noisy data. *Inverse Problems*, 26(9):095016, 2010.

[11] Laurent Bourgeois and Jérémi Dardé. The morozov's principle applied to data assimilation problems. *ESAIM: M2AN*, 56(6):2021-2050, 2022.

[12] Laurent Bourgeois and Arnaud Recoquillay. A mixed formulation of the Tikhonov regularization and its application to inverse PDE problems. *ESAIM Math. Model. Numer. Anal.*, 52(1):123-145, 2018.

[13] Haïm Brezis. *Analyse fonctionnelle : théorie et applications*. Editions Dunod, 1999.

[14] Pablo J. Cobelli, Philippe Petitjeans, Agnès Maurel, and Vincent Pagneux. Determination of the bottom deformation from space- and time-resolved water wave measurements. *Journal of Fluid Mechanics*, 835:301-326, 2018.

[15] David Colton and Rainer Kress. *Inverse acoustic and electromagnetic scattering theory*, volume 93. Springer, 1998.

[16] Pierre Demeestère. A remark on the relation between the tykhonov regularization and constraint relaxation for an optimal control problem. *Applied Mathematics Letters*, 11(5):85-89, 1998.
[17] Juliette Dubois, Sébastien Imperiale, Anne Mangeney, François Bouchut, and Jacques Sainte-Marie. Acoustic and gravity waves in the ocean: a new derivation of a linear model from the compressible euler equation. *Journal of Fluid Mechanics*, 970:A28, 2023.

[18] Denis Dutykh. *Modélisation mathématique des tsunamis*. PhD thesis, Ecole Normale Supérieure de Cachan, 2007.

[19] Ivar Ekeland and Roger Temam. *Convex analysis and variational problems*. SIAM, 1999.

[20] Alexandre Ern and Jean-Luc Guermond. *Theory and practice of finite elements*, volume 159. Springer, 2004.

[21] Lawrence Craig Evans and Ronald Francis Gariepy. *Measure theory and fine properties of functions*. Studies in advanced mathematics. CRC Press Boca Raton, FL, 1992.

[22] Marco A. Fontelos, Rodrigo Lecaros, Juan López-Ríos, and Jaime H. Ortega. Bottom detection through surface measurements on water waves. *SIAM Journal on Control and Optimization*, 55(6):3890–3907, 2017.

[23] David Gilbarg and Neil Trudinger. *Elliptic Partial Differential Equations of Second Order*. Grundlehren der mathematischen Wissenschaften. Springer Berlin Heidelberg, 2013.

[24] Bernabe Gomez and Usama Kadri. Earthquake source characterization by machine learning algorithms applied to acoustic signals. *Scientific Reports*, 11(1):23062, 2021.

[25] Risto Heijmans. When does the expectation of a ratio equal the ratio of expectations? *Statistical Papers*, 40(1):107–115, 1999.

[26] Sergey Kabanikhin and Olga Krivorotko. A combined inverse tsunami problem. *Bull. Nov. Comp. Center, Math. Model. in Geoph*, 16:45–58, 2013.

[27] N. Kielbasiewicz and E. Lunéville. User documentation for xlife++. https://uma.ensta-paris.fr/soft/XLiFE++/var/files/docs/usr/user documentation.pdf, 2019.

[28] Andreas Kirsch. *An Introduction to the Mathematical Theory of Inverse Problems*. Applied Mathematical Sciences. Springer International Publishing, 2021.

[29] David Lannes. *The water waves problem: mathematical analysis and asymptotics*, volume 188. American Mathematical Soc., 2013.
[30] Marc Lenoir and A. Tounsi. The localized finite element method and its application to the two-dimensional sea-keeping problem. *SIAM Journal on Numerical Analysis*, 25(4):729-752, 1988.

[31] Jacques-Louis Lions and Enrico Magenes. *Problèmes aux limites non homogènes et applications. Vol. 2*. Travaux et Recherches Mathématiques, No. 18. Dunod, Paris, 1968.

[32] Lecaros Rodrigo, Juan López-Ríos, Jaime H. Ortega, and Sebastian Zamorano. The stability for an inverse problem of bottom recovering in water-waves. *Inverse Problems*, 36(11):115002, oct 2020.

[33] Jerome A. Smith. Revisiting oceanic acoustic gravity surface waves. *Journal of Physical Oceanography*, 45(12):2953–2958, 2015.

[34] Michael Stiassnie. Tsunamis and acoustic-gravity waves from underwater earthquakes. *Journal of Engineering Mathematics*, 67:23-32, 2010.

[35] Vishal Vasan and Bernard Deconinck. The inverse water wave problem of bathymetry detection. *Journal of Fluid Mechanics*, 714:562-590, 2013.

[36] Tokuo Yamamoto. Gravity waves and acoustic waves generated by sub-sea earthquakes. *International Journal of Soil Dynamics and Earthquake Engineering*, 1(2):75-82, 1982.